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Preface 


The  purpose  of  this  thesis  was  to  examine  electron  kinetics  in  plasmas.  The  conditions 
encompass  those  relevant  to  the  development  of  both  excimer  lasers  and  negative  ion  sources. 
The  vehicle  for  this  investigation  was  a  numeric  solution  of  the  time-dependent  collisional  Boltz¬ 
mann  equation.  The  method  of  solution  is  based  upon  that  developed  by  Rockwood  (14).  The 
processes  considered  in  this  code  are  elastic  electron-neutral  collisions,  as  well  as  ionization, 
excitation,  attachment,  and  superelastic  collisions. 

Although  much  re-structuring  and  many  program  modifications  and  additions  were  required 
in  order  to  study  these  devices  and  include  these  processes,  the  code’s  foundation  was  laid  by 
Ges  Seger.  His  programming  style  and  foresight  made  it  much  easier  for  me  to  pick  up  where  he 
left  off.  I  also  wish  to  thank  my  advisor,  Dr  William  Bailey  for  his  insight,  help,  and  suggestions 
throughout  this  project.  He  even  went  as  far  as  France  in  order  to  ask  my  questions  of  Dr  Bre¬ 
tagne,  another  person  that  I  am  indebted  to.  His  many  articles  on  the  subjects  pertaining  to 
plasma  physics  gave  me  further  insight  and  ideas  that  proved  useful. 

Finally,  I  would  like  to  thank  my  wife  Kathi  for  her  patience  and  support  throughout  this 
project,  as  well  as  staying  up  with  me  to  the  wee  hours  needlepointing  beside  my  computer,  just 
so  we  could  be  together. 
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Abstract 


In  a  continuing  effort  to  analyze  plasmas  generated  in  a  variety  of  ways,  several  refinements 
and  improvements  have  been  added  to  a  numerical  solution  of  the  time-dependent  collisional 
Boltzmann  equation.  The  computational  method  utilized,  originated  by  Seger  (16)  with  the 
Rockwood  formalism  (14),  includes  elastic  collisions,  excitation,  ionization,  and  superelastic  col¬ 
lisions  from  multiple  states,  as  well  as  electron  attachment.  The  attachment  process  is  explored, 
with  the  results  confirming  predictions  made  by  the  energy  moment  of  the  Boltzmann  equation. 
The  possibility  of  steady  state  negative  mobility  is  explored  in  lean  Ar/F2  mixes,  resulting  in 
computational  confirmation  of  theoretical  predictions.  A  non-uniform  energy  grid  is  also  incor¬ 
porated  into  the  method  of  solution,  requiring  new  finite  differenced  representations  for  the  elec¬ 
tron  flux  terms  from  field,  recoil,  and  electron-electron  interactions.  All  of  these  terms  are  tested 
independently  and  validated  against  analytic  results.  Implementation,  complete  with  inelastic 
collisions,  is  validated  against  experimental  transport  data  for  both  molecular  and  rare  gases.  A 
method  for  dynamically  allocating  the  energy  grid  in  an  effort  to  optimize  the  computation  is 
developed  and  evaluated.  The  inclusion  of  the  non-uniform  grid  results  in  the  consideration  of 
larger  energy  ranges  than  previously  possible.  This  larger  energy  range,  demonstrated  up  to  90 
eV,  allows  the  exploration  of  the  effects  of  electron  beams  on  the  electron  energy  distribution. 
Energy  distribution  functions  typical  of  a  10  A  electron  beam  source  are  calculated  at  60  and  90 
eV  and  compared  to  previous  solutions. 


ELECTRON  KINETIC  STUDIES  UTILIZING  THE  SOLUTION  OF  THE 
TIME-DEPENDENT  COLLISIONAL  BOLTZMANN  EQUATION 


I.  Introduction 

The  Boltzmann  equation  is  a  continuity  equation  describing  the  flow  of  particles  through 
phase  space  under  the  influence  of  externally  applied  forces  and  collisions.  While  analytic  solu¬ 
tions  to  the  Boltzmann  equation  exist  for  special  casus,  a  quantitative  analysis  of  plasma  devices 
demands  that  the  solution  be  obtained  through  the  use  of  numerical  methods.  While  some 
methods  rely  on  finding  only  the  steady  state  solution  (8,10),  a  more  general  method  is  adapted. 
This  admits  the  study  of  devices  such  as  pulsed  electron  beam  sustained  plasmas,  and  permits 
self-consistent  coupling  of  plasma  chemistry  and  heavy  particle  kinetics.  This  time-dependent 
solution  can  be  propagated  forward  in  time,  allowing  the  observation  of  the  temporal  evolution 
of  the  distribution,  until  either  a  steady  state  solution  is  achieved  or  a  self-similar  form  of  the 
distribution  is  obtained. 

Solutions  of  the  Boltzmann  equation  in  energy  space  yield  electron  Energy  Distribution  Func¬ 
tions  (EEDF).  Once  the  EEDF  is  obtained,  transport  parameters  and  plasma  kinetic  rates  can 
then  be  calculated  by  taking  various  moments  of  the  distribution.  These  parameters  and  rates  are 
invaluable  in  the  design  and  optimization  of  many  plasma  devices.  In  gas-discharge  lasers,  a 
knowledge  of  the  distribution  function  allows  the  kinetic  pumping  rates  to  be  determined,  which 
are  a  critical  part  of  laser  design  and  optimization.  Since  the  plasma  within  such  a  device  is  gen- 
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erated  and  sustained  by  electron  interactions  with  the  electric  field  and  gas  particles,  accurate 
pumping  calculations  are  critically  dependent  on  the  proper  treatment  of  all  relevant  collisional 
processes.  Such  processes  include  impact  ionization,  excitation,  and  superelastic  collisions. 

In  excimer  lasers  and  generators,  ancilliary  sources  of  ionization  in  the  form  of  an  electron 
beam  or  hot  filament  are  employed.  Beam  energies  in  such  devices  are  typically  at  much  higher 
energies  than  the  plasma  mean  energy.  Therefore,  the  analysis  of  these  devices  requires  not  only 
the  inclusion  of  a  source  term,  but  also  an  extension  of  the  energy  range  in  the  solution  method. 
In  the  exploration  of  the  excimer  lasers,  and  in  particular  the  e-  beam  sustained  rare  gas  halide 
lasers,  the  process  of  electron  attachment  is  of  fundamental  importance,  since  it  is  just  such  a 
process  that  allows  the  production  of  the  negative  ions  that  are  used  in  the  reaction  mechanism. 
Thus,  any  analysis  of  EEDFs  characteristic  of  excimer  lasers  must  properly  include  the  process 
of  electron  attachment.  The  process  of  attachment  is  also  of  fundamental  importance  in  dis¬ 
charges  that  contain  a  halogen  gas  as  one  of  the  constituents.  These  gases  are  very  electronega¬ 
tive,  and  as  such  will  tend  to  remove  free  electrons  from  the  system.  This  electron  loss  can  result 
in  unusual  effects  within  the  discharge.  One  such  effect  leads  to  the  possibility  of  a  plasma 
discharge  with  a  steady  state  negative  total  elecmjn  mobility.  When  an  electric  field  is  applied  in 
a  typical  discharge,  the  electrons  will  tend  to  drift  in  a  direction  anti-parallel  to  the  field.  In  such 
a  case,  the  mobility  is  said  to  be  positive.  A  negative  mobility,  therefore,  leads  to  the  electrons 
drifting  parallel  to  the  electric  field.  In  such  a  case,  it  is  possible  for  the  electrons  to  give  energy 
to  the  field,  resulting  in  an  increase  in  the  field  strength  (15:1596). 

The  purpose  of  this  study  was  to  explore  the  electron  processes  fundamental  to  the  operation 
of  plasma  devices  that  are  of  interest  to  the  Air  Force.  These  devices  include  electron-beam  sus¬ 
tained  plasmas,  excimer  lasers,  and  magnetic  multicusp  discharges.  The  analysis  of  each  of 
these  requires  additional  consideration  within  the  code  in  terms  of  energy  conservation,  external 
sources,  and  electron  interactions,  both  between  themselves  and  with  the  gas  particles.  All  of 
these  processes  will  contribute  to  the  steady  state  operating  point  of  the  plasma,  and  thereby  have 
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a  very  large  impact  on  the  design,  optimization,  and  analysis  of  the  devices. 

The  analysis  and  understanding  of  the  physics  associated  with  these  devices  can  be  greatly 
improved  by  careful  calculation  of  those  properties  which  can  be  experimentally  measured  in  the 
laboratory.  If  these  calculated  parameters  compare  favorably  to  the  experimentally  measured 
ones,  then  confidence  can  be  given  that  all  pertinent  electron  processes  have  been  identified  in 
the  solution  method.  Such  parameters  include  drift  velocity,  average  energy,  characteristic 
energy,  and  the  ratio  of  the  Townsend  alpha  coefficient  to  the  neutral  number  density,  a/N. 

In  the  analysis  of  excimer  lasers,  multicusp  discharges,  and  electron  beams,  a  source  of  elec¬ 
trons  is  typically  present  at  energies  which  are  high  compared  to  the  energies  typical  of  gas  dis¬ 
charges.  .Therefore,  in  order  to  study  such  devices,  the  energy  range  of  the  solution  must  be 
extended  to  energies  in  the  range  of  0.1-100  keV.  In  the  method  of  solution  developed  by 
Rockwood  (14),  this  energy  range  is  broken  into  K  energy  bins  of  equal  size,  and  the  method  of 
finite  differences  utilized  to  solve  the  resulting  K  coupled  equations  simultaneously.  This  uni¬ 
form  energy  method  is  practical  in  devices  in  which  the  energy  range  can  be  broken  into  approxi¬ 
mately  300  energy  bins  or  less.  In  the  solution  method,  the  K  coupled  equations  are  arranged 
into  a  square  matrix  of  size  K  x  K,  and  the  inverse  calculated.  Aside  from  the  physical  limitation 
on  the  size  of  the  K  matrix  due  to  computer  memory,  the  time  required  for  this  inverse  to  be 
calculated  is  proportional  to  the  cube  of  K  (12:30-31, 33,  368),  regardless  of  the  method  used. 
This  requirement  translates  into  a  maximum  energy  range  of  approximately  30  eV  for  rare  gases, 
and  15  eV  or  so  for  molecular  gases.  The  number  of  energy  bins  required  to  accurately  represent 
the  distribution  without  rendering  the  basic  assumption  of  finite  differencing  invalid  is  critically 
linked  to  the  cross  section’s  dependence  on  energy.  In  molecular  gases,  the  inelastic  processes 
that  are  grouped  near  the  lower  energies  require  a  small  energy  bin  size  in  order  to  accurately 
represent  the  distribution.  This  smaller  bin  size  reduces  the  energy  range  of  the  solution  for  a 
fixed  number  of  bins.  In  rare  gases  the  threshold  energy  of  the  inelastic  processes  are  higher, 
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and  a  larger  bin  size  can  generally  be  used  before  it  leads  to  invalidating  the  finite  differenced 
assumption.  Clearly,  then,  this  uniform  bin  size  cannot  be  efficiently  used  to  extend  the  energy 
range  to  that  required  by  the  e-  beam  devices. 

In  order  to  extend  the  present  method  of  solution  of  the  Boltzmann  equation  out  to  these 
higher  energies,  while  still  retaining  accuracy  and  stability  of  the  solution,  several  researchers 
have  adopted  a  non-uniform  energy  grid,  consisting  of  differential  energy  elements  of  unequal 
size  (1:1201)  (2:2209).  One  method  used  is  to  allocate  the  energy  axis  logarithmically  (3:814). 
This  method  relies  on  the  premise  that  the  distribution  changes  most  rapidly  at  the  lower  ener¬ 
gies,  and  less  so  as  the  energy  is  increased.  This  method  can  be  effective  at  calculating  distrib¬ 
utions  out  to  the  higher  energies  as  Bretagne,  et  al  (3)  has  done  for  molecular  hydrogen,  in  which 
the  inelastic  thresholds  are  largely  at  the  lower  energies  (1/2  to  12  eV).  However,  this  method 
does  not  favor  gases  with  inelastic  processes  which  occur  at  higher  energies,  such  as  for  the  rare 
gases.  This  is  because  the  differential  energy  elements  are  always  increasing  with  energy,  with 
no  regard  as  to  where  an  inelastic  process  begins.  Since  inelastic  processes  can  lead  to  signifi¬ 
cant  structure  in  the  distribution,  this  method  is  not  as  desirable  as  one  in  which  the  energy  ele¬ 
ment  size  is  reduced  in  regions  around  the  threshold  energy. 

While  the  non-uniform  energy  bin  method  is  used  to  extend  the  energy  range  of  the  calcula¬ 
tion,  it  may  also  prove  to  be  more  accurate  than  the  uniform  grid  in  the  computation  of  the  solu¬ 
tion  for  a  fixed  energy  range.  Thus,  if  there  is  some  method  which  can  increase  the  accuracy  of 
the  solution  in  the  same  amount  of  time,  then  clearly  it  is  advantageous  to  include  this  method 
into  the  code. 

A  computer  code  that  yields  the  numerical  solution  to  the  time-dependent  Boltzmann  equa¬ 
tion  along  a  non-uniform  energy  grid  has  been  developed  using  the  formalism  given  by  Rock- 
wood  (14),  with  the  foundation  established  by  Seger  (16).  This  code  considers  excitation, 
ionization,  attachment,  elastic  and  superelastic  processes  as  well  as  electron-electron  interactions 
and  momentum  transfer.  The  possibility  of  excitation,  ionization,  and  superelastic  collisions 
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from  multiple  excited  states  is  included.  The  transport  parameters  calculated  with  the  distrib¬ 
utions  resulting  from  the  non-uniform  finite  differenced  representations  of  the  field  and  recoil 
driven  flux  terms  are  validated  against  analytic  results.  The  new  finite  differenced  terms 
representing  electron  -  electron  collisions  are  tested  for  conservation  of  energy  and  particles,  as 
well  as  for  convergence  to  the  proper  form  when  considered  alone.  Using  a  logarithmically  allo¬ 
cated  energy  grid,  transport  parameters  for  Ar  and  H2  are  compared  to  uniform  grid  calculations 
of  the  same,  as  well  as  experimental  data.  EEDFs  for  an  electron  beam  generated  H2  plasma  are 
calculated  using  a  logarithmic  energy  grid  and  compared  to  previous  calculations.  A  new 
method  has  been  developed  which  allocates  the  energy  axis  based  on  the  energy  dependence  of 
elastic  and  inelastic  processes.  The  utility  of  this  dynamic  energy  bin  allocation  method  is 
explored  using  an  approximate  solution  to  the  time  independent  Boltzmann  equation  as  a  bench¬ 
mark  distribution.  In  addition,  the  process  of  electron  attachment  is  tested  and  validated  with 
predictions  made  by  the  energy  moment  of  the  Boltzmann  equation.  The  possibility  of  steady 
state  negative  mobility  was  studied  in  lean  Ar/F2  gas  mixtures,  with  the  results  confirming  theo¬ 
retical  predictions. 
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Origins  of  Solutions  to  the  Boltzmann  Equation 


The  methods  used  to  solve  the  Boltzmann  equation  for  the  distribution  of  electrons  in  energy 
space  date  back  to  1946,  with  T.  Holstein  (8).  Although  his  equations  were  derived  for  distrib¬ 
utions  driven  by  AC  field  conditions,  they  can  be  applied  to  the  DC  case  as  well.  Holstein 
started  with  the  Boltzmann  equation  describing  the  continuity  of  electrons  in  phase  space: 


|-(v.Vl)f+(?-Vv)f=(f 


(1) 


/collisions 


where  f  is  the  electron  distribution  function  describing  the  number  of  electrons  in  six  dimen¬ 
sional  phase  space,  thus  /= /(vJf,v>,v„x,y,z,r),  v  is  the  electron  velocity  vector,  a  is  the  elec¬ 
tron  acceleration  due  to  some  applied  force,  V,  is  the  gradient  in  geometric  coordinates,  and  Vv  is 
the  gradient  in  velocity  coordinates.  Holstein  then  assumed  that  the  function  f  is  very  nearly 
spherically  symmetric  in  velocity  space  for  a  sufficiently  small  applied  field.  This  assumption 
allows  the  expansion  of  the  distribution  into  an  infinite  series  of  spherical  harmonics,  with  the 
axis  defined  by  the  direction  of  the  field.  Assuming  azimuthal  symmetry,  the  spherical  harmon¬ 
ics  move  to  Legendre  polynomials: 


f(x,v,Q,t)  =  I  /„(x,v,OP„(cose) 

n  =0 


(2) 


where  the  P„’s  are  the  Legendre  polynomials.  In  practice,  this  infinite  series  is  truncated  after 

the  first  two  terms,  and  is  called  the  Lorentz,  or  P„  approximation.  The  result  of  this  expansion 
is  a  pair  of  coupled  equations  in  terms  of  f0  and/,,  where  the  former  represents  the  symmetric 


part  of  the  distribution,  while  the  latter  represents  the  anisotropic  part.  At  this  point,  Holstein 
grouped  the  pair  of  coupled  equations  into  a  single  equation.  He  then  concentrated  on  finding 
the  steady  state  solution  to  this  single  equation,  which  represents  the  solution  to  the  Boltzmann 
equation  for  D.C.  or  high  frequency  A.C.  discharges. 

Computational  Method  of  Solution 

In  1973,  Stephen  Rockwood  published  a  paper  outlining  a  method  used  to  calculate  time 
dependent  numerical  solutions  (14).  In  this  method,  the  processes  considered  were  electron- 
neutral  elastic  interactions,  inelastic  interactions  in  the  form  of  excitation,  ionization,  and  super¬ 
elastic  collisions,  as  well  as  electron-electron  collisions.  The  basic  principle  used  to  describe  the 
time  rate  of  change  of  the  electron  number  density  was  the  continuity  equation.  This  equation 
can  be  written  as 

^+V-J  =  S-L  (3) 

at 

where  n  is  the  number  density  of  electrons,  and  V »  J  represents  the  flux  of  those  electrons 

through  energy  space  due  to  elastic  and  inelastic  collisions  as  well  as  the  external  force  arising 
from  the  electric  field.  S  and  L  represent  an  external  source  and  loss  term  respectively.  Thus, 
the  time  rate  of  change  of  the  electron  number  density  due  to  elastic  and  inelastic  collisions,  as 
well  as  external  forces  can  be  written  as  a  sum  of  terms  expressing  the  divergence  of  the  flux  of 
electrons  through  energy  space  resulting  from  each  process: 
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+  S-L 


(4) 


dn  _  dJ/  <Hei  diee  ^  f 

3 dt  +  U' )collisions 

Here,  the  subscripts  f,  el,  and  ee  denote  the  field,  elastic,  and  electron  -  electron  terms  respec¬ 
tively.  As  Rockwood  shows  (14:2348),  the  current  density  of  electrons  in  energy  space  due  to 
the  applied  external  electric  field  can  be  written  as 

2Ne2(E/N)2ef  n__ dn' 
f  3m(v/N)  ^2e  de  } 

v 

-=  -  I  gsa,(e)  (5) 

N  \  m  J  s 

where  v  and  a(e)  are  the  collision  frequency  and  cross  section  for  momentum  transfer  respec¬ 
tively,  q,  is  the  mole  fraction  of  species  s,  and  n  is  the  electron  number  density  at  energy  e. 

Similarly,  Rockwood  shows  that  the  current  density  of  electrons  in  energy  space  due  to  elastic 
collisions  with  molecules  of  species  s  at  a  gas  temperature  of  T  is  (14:2348) 

-f  f  kT  ^  3rt" 

Jel=  L  6 

V  =  2mN(2EJm)mlqsos(e)/Ms 

I 

where  M,  is  the  molecular  weight  of  species  s,  and  k  is  Boltzmann’s  constant.  Note  that  in  both 

of  the  expressions  for  the  current  densities  J j  and  J,,  there  is  a  term  that  is  proportional  to  the 
number  density  of  electrons  n  (convective),  and  another  which  is  proportional  to  the  energy  gra¬ 
dient  of  the  number  density  dn/de,  (diffusive). 
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The  remaining  terms  in  the  continuity  (Boltzmann)  equation  are  those  which  represent  inelas¬ 
tic  electron  interactions  with  gas  molecules  and  electron  -  electron  interactions.  The  inelastic 
terms  can  be  written  classically  as  n,No(e)v(e),  where  the  cross  section  a  would  be  for  the 
appropriate  inelastic  process,  N  represents  the  number  of  collision  partners,  while  n,  represents 
the  number  of  electrons  at  the  collision  energy. 

Intra-electron  interactions  can  also  be  represented  by  the  use  of  the  Fokker-Planck  approxi¬ 
mation  to  the  electron-electron  collision  integral  described  by  the  energy  space  divergence  of  the 
electron  current  density  (14:2357): 


3*, 

3/(E,,)  3e 


J„  =  al 


r  ( 
P 
L  V 


n  dn 

2e  de 


\ 


-Qn 


where 


(7) 


P(e,0=2e  1/2  Jxn(x,t)dx 

o 

oo 

+  2e  Jx~]/2n(x,t)dx 

t 

e 

Q(e,f)  =  J n(x,t)dx 
o 

Following  the  traditional  method  of  solving  differential  equations  by  the  method  of  finite 
differencing,  Eq  (4)  is  projected  onto  an  energy  axis  that  has  been  divided  into  K  energy  bins  of 
in  general,  various  widths.  Thus  the  differential  equation  is  recast  into  K  coupled  linear  equa- 
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tions  and  solved  simultaneously.  The  actual  derivation  of  the  finite  differenced  terms  represent¬ 
ing  the  flux  of  electrons  through  energy  space  due  to  the  applied  electric  field  and  elastic 
collisions  between  electrons  and  neutral  gas  molecules  is  included  in  Appendix  A.  In  this 
derivation,  particular  attention  has  been  paid  to  the  fact  that  the  energy  axis  is  constructed  of 
energy  bins  of  various  sizes,  as  well  as  maintaining  numerical  stability  with  no  electric  field 
present.  Neglecting  electron  -  electron  collisions  due  to  their  non-linearity,  and  setting  to  zero 
the  source  and  loss  terms,  the  finite  differenced  Boltzmann  equation  can  be  written  as  (14:2349) 


dn 


—  =  ak  _  }nk  _  i  +  bk  +  }nk  + 1  -  (ak  +  bk)nk  +  I  N,(R^ +mnk +m  .  +  Rsjk_mNJJNs 


S'J 


+Rlt  +msink  +m„  +  SU  ^  Kmnm  ~  (R#  +  Rsjk  +  RltK) 


(B) 


where  k  is  the  energy  bin,  s  is  the  gas  species,  j  identifies  a  particular  inelastic  process,  and  m  is 
the  threshold  energy  for  that  process.  The  R’s  represent  rates  (a(e)v(e))  for  each  process.  The 
ak s  and  bk' s  represent  elastic  promotions  and  demotions  respectively,  each  from  the  k* 
energy  bin  due  to  the  applied  field  and  recoil  with  heavy  particles,  while  nk  represents  the  num¬ 
ber  density  of  electrons  with  energy  between  et_!  and  £*.  Therefore,  nk  is  a  centered  quantity  in 
relation  to  the  bin  spacing.  Thus  the  first  term  on  the  right  represents  the  upflux  of  electrons  due 
to  momentum  transfer  and  the  external  field  from  the  energy  bin  directly  below  the  k*  bin.  Simi¬ 
larly,  the  second  term  represents  the  downflux  of  electrons  due  to  momentum  transfer  and  field 
from  the  energy  bin  directly  above  the  k^bin.  The  third  term  represents  a  loss  of  electrons  out  of 
the  k*  bin  due  to  the  upflux  of  electrons  to  bin  k+1,  and  downflux  of  electrons  to  bin  k-1.  The 
first  term  in  the  summation  represents  an  addition  of  electrons  to  bin  k  due  to  an  electron  that 
underwent  impact  excitation  at  an  energy  level  of  k+m,  where  m  represents  the  threshold  energy 
for  the  excitation  process.  When  that  electron  struck  an  atom  and  excited  it  to  an  upper  level,  it 
expended  an  energy  equal  to  the  excitation  energy  for  that  process.  Thus,  the  electron  lost  an 
amount  of  energy  equal  to  the  threshold  energy,  thereby  becoming  an  addition  to  the  energy  bin 
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that  is  located  a  "distance"  of  the  threshold  energy  below  it  on  the  energy  axis.  The  second  term 
in  the  summation  represents  superelastic  contributions  to  the  k*  bin  from  the  energy  bin  located  a 
distance  of  the  threshold  energy  below  it.  Likewise,  the  third  term  represents  additions  to  bin  k 
due  to  electrons  causing  impact  ionization  at  bin  k+m,  where  m  represents  the  threshold  for  the 
ionization  process.  The  next  term  reflects  the  assumption  that  all  of  the  secondary  electrons 
emitted  by  ionization  reappear  in  the  first  energy  bin  on  the  axis.  The  last  terms  under  the  sum¬ 
mation  represent  all  the  inelastic  losses  out  of  the  k*  bin  -  excitation,  superelastic,  and  ionization, 
respectively. 

Recognizing  that  this  differential  equation  was  constructed  so  as  to  be  linear  in  nk,  it  can  be 
rewritten  in  matrix  form  as 


drt* 

d / 


~  Q.  lnl 


(9) 


where  C  is  a  K  x  K  matrix  that  remains  constant  with  time,  with  K  denoting  the  number  of 
energy  bins  along  the  entire  energy  axis.  The  applied  electric  field  and  recoil  terms  will  be 
placed  into  C  such  that  it  will  become  tridiagonal.  The  excitation  and  ionization  terms  will  be 
placed  into  the  C  matrix  above  the  diagonal,  while  the  elements  representing  superelastic  con¬ 
tributions  will  be  placed  below  the  diagonal.  Thus,  every  element  along  the  diagonal  itself  repre¬ 
sents  a  loss  out  of  the  k*  bin,  while  every  element  not  on  the  diagonal  will  represent  an  influx  of 
electrons  into  the  k*  bin. 

The  left  hand  side  of  Eq  (9)  can  be  differenced  using  a  simple  Euler  relation: 


dnk  nk(t  + At)- nk(t) 
d \t~  At 


(10) 
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Since  the  n,  on  the  right  hand  side  of  Eq  (9)  is  the  electron  density  at  time  t  +  Ai,  Eq  (9)  can  be 


rewritten  as 


n(t+At)  =  (I-CA/)_1«(0  (11) 

where  I  is  the  identity  matrix,  and  7i(t)  represents  the  electron  number  density  in  each  of  the  K 
energy  bins  at  time  t.  From  Eq  (1 1),  it  can  be  seen  that  if  the  electron  number  density  is  known 
at  any  time  t,  it  can  easily  be  calculated  for  the  time  t  +  At.  Since  the  C  matrix  is  a  constant,  the 
inverse  need  be  computed  only  once.  Rockwood  states  that  the  convergence  properties  of  this 
solution  method  are  very  good,  even  for  arbitrary  choices  of  n*(0)  (14:2349). 

With  the  C  matrix  incorporating  all  of  the  electron-neutral  collisional  interactions,  the  single 
differential  equation  Eq  (1)  is  broken  into  a  number  of  coupled  equations.  This  system  of  equa¬ 
tions  is  then  solved  simultaneously,  with  the  solution  being  given  by  the  distribution  function  at 
the  next  time  step.  This  solution  is  then  marched  forward  in  time,  until  the  distribution  changes 
very  little  in  comparison  to  the  previous  solution.  At  such  a  point  the  self-similar  form  of  the 
distribution  has  been  achieved.  This  solution  method  has  not,  however,  included  the  effects  of 
electron-electron  collisions  in  determining  the  final  distribution.  In  the  absence  of  any  other  pro¬ 
cesses,  the  electron-electron  collisions  will  drive  the  distribution  to  a  Maxwellian.  The  influence 
of  this  effect  becomes  more  important  at  the  higher  fractional  ionizations.  However,  since  they 
are  non-linear  they  cannot  be  included  into  the  C  matrix,  but  must  be  handled  another  way. 

Electron-Electron  Collisions 

The  electron-electron  interactions  are  nonlinear  with  respect  to  n,  thus  they  cannot  be 
included  into  the  C  matrix  of  Eq  (9).  They  are  incorporated  into  the  solution  method  by  adding 
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to  the  right  hand  side  of  Eq  (9)  the  matrix  T(n ),  which  is  tridiagonal.  The  resulting  equation  can 
be  written  as: 


n(t+Ai)  =  (I-CA/)'1  (I  +  TA/)7f(r)  (12) 

The  elements  within  T  directly  correspond  to  the  C  matrix  for  the  special  case  of  momentum 
transfer  only.  The  elements  along  the  diagonal  represent  losses  out  of  the  bin,  while  the  elements 
above  and  below  the  diagonal  represent  the  downflux  and  upflux,  respectively,  of  electrons  from 
other  bins.  The  primary  difference  between  the  two  treatments  lies  in  the  fact  that  the  coulomb 
collision  integral  takes  the  entire  distribution  of  electrons  into  account  at  each  time  step  in  order 
to  describe  the  electron  interactions,  which  is  the  source  of  the  nonlinearity.  The  momentum 
transfer  between  electrons  and  neutral  particles,  on  the  other  hand,  requires  only  a  simple  binary 
collision,  thus  they  can  be  treated  implicitly,  while  the  former  is  treated  explicitly  in  this 
approach.  The  elements  of  the  T  matrix  must  be  computed  for  each  time  iteration  throughout  the 
computation.  Appendix  B  contains  the  derivation  of  the  finite  differenced  terms  representing 
electron -electron  interactions  along  the  non-uniform  energy  axis.  In  this  derivation,  it  is  impor¬ 
tant  to  remember  that  the  finite-differenced  representation  for  the  electron-electron  interactions 
must  obey  three  important  principles:  conservation  of  energy,  conservation  of  particles,  and 
relaxation  of  the  distribution  toward  a  Maxwellian  when  electron-electron  collisions  are  the  only 
process  at  work.  Conservation  of  particles  is  ensured  by  the  flux  divergent  form  of  the  bin  rate 
equations,  coupled  with  appropriate  boundary  conditions.  Conservation  of  energy  is  guaranteed 
by  proper  construction  of  the  electron-electron  matrices.  The  last  property  requires  that  with 
electron-electron  interactions  only,  the  distribution  of  particles  should  become  a  Maxwellian  at  a 
temperature  equal  to  2/3  of  the  average  energy.  This  property  is  achieved  by  the  use  of  detailed 
balance  in  requiring  that  the  upflux  of  particles  from  one  bin  to  the  next  higher  bin  equals  the 
downflux  from  the  higher  bin  to  the  one  just  below  it. 
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Seger  -  A  Working  Code 


Developed  using  the  Rockwood  formalism,  Seger’ s  code  (16)  was  able  to  solve  the  Boltz¬ 
mann  equation  for  the  self-similar  form  of  the  solution.  In  this  method,  the  n(t  +  At)  calculated 
by  Eq  (11)  becomes  the  next  7T(r)in  the  process.  This  code  considered  momentum  transfer, 
superelastic,  excitation,  and  ground  state  ionization  processes.  What  was  required  of  the  user 
were  input  tables  of  cross  section  data  for  each  process  of  each  gas  considered,  the  relative 
amount  of  each  gas  present,  a  particular  E/N  value,  and  certain  other  parameter  which  defined 
the  problem  to  solve,  such  as  neutral  gas  temperature  and  pressure,  excited  state  population,  etc. 
Once  the  solution  to  the  Boltzmann  equation  was  obtained,  then  drift  velocity,  average  energy, 
and  the  free  diffusion  coefficient  could  be  calculated. 


Transport  Parameters 


The  transport  parameters  serve  to  link  L.c  microrccnic  world  to  the  macroscopic.  In  practice 
it  is  very  difficult  for  the  experimentalist  to  measure  directly  the  distribution  function  itself,  but  it 
is  very  practical  to  measure  other  parameters  which  in  theory  depend  directly  on  the  distribution 
function.  Conversely,  the  theoretician  cannot  directly  determine  any  of  the  transport  parameters, 
but  must  first  calculate  the  distribution  function  upon  which  they  are  based.  In  order  to  deter¬ 
mine  if  this  distribution  function  is  valid,  the  parameters  calculated  from  it  must  be  compared  to 
the  experimentally  measured  ones.  In  a  very  real  sense  then,  the  experimentalist  is  much  like  a 
blind  man  feeling  different  parts  of  the  elephant  in  an  effort  to  picture  the  whole  beast. 

The  calculation  of  the  transport  parameters  is  based  upon  taking  different  moments  of  distri¬ 
bution  function.  The  drift  velocity  is  calculated  by  (10:472): 
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(13) 


o 

where  the  two  term  approximation  has  been  explicity  retained  in  this  expression.  Thus, 
f(u)=  f0(u)  +  /, (u) cos 0  ,  where /e  and  /,  are  the  reduced  distributions  and  are  coupled  to 
each  other.  The  free  diffusion  coefficient  D/  is  calculated  by  (17:176): 


D  = 


■n 

3^ 


(14) 


where  u  is  the  velocity,  vr  is  the  system  relaxation  collision  frequency  and  is  equal  to  the  effec¬ 
tive  collision  frequency  when  v  is  independent  of  velocity  (17:171-176).  The  average  energy  is 
calculated  by  taking  the  energy  moment  of  the  solution  to  the  Boltzmann  equation: 


oo 

J 


e/(e)de 


<£>  = 


J- 


(15) 


/(e)  de 


The  characteristic  energy  is  a  quantity  which  can  be  measured  in  an  actual  plasma  discharge, 
and  it  can  be  used  to  define  an  equivalent  Maxwellian  temperature  (4:101): 


£*  = 


eD 


=  kT 


(16) 
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Here,  it  has  been  written  in  terms  of  the  Einstein  relationship.  It  is  seen  that  if  the  distribution 
was  indeed  Maxwellian,  then  the  characteristic  energy  would  be  2/3  of  the  average  energy. 

Moments  of  the  Boltzmann  Equation 

Serving  as  another  link  between  the  unseen  world  of  the  very  small  and  the  world  in  which 
we  live  are  the  moments  of  the  Boltzmann  equation.  Describing  different  aspects  of  the  plasma, 
the  moments  of  the  Boltzmann  equation  detail  the  macroscopic  quantities  that  researchers  are 
able  to  measure.  These  are,  principally,  the  production  and  loss  of  electrons,  the  velocity  of  elec¬ 
trons,  and  the  energy  of  the  electrons.  The  moments  of  the  Boltzmann  equation  are  generally 
written  in  terms  of  the  averages  which  emphasize  their  macroscopic  viewpoint. 

The  continuity  equation  has  been  mentioned  previously  in  some  detail,  but  will  be  written 
here  in  terms  of  the  coordinate  space  divergence  of  the  electron  flux: 

—  =  V  •  nv  +  S  -  L  (17) 

d t 

Assuming  that  the  plasma  is  spatially  homogeneous,  and  rewriting  the  production  and  loss  terms 
as  frequencies,  Eq  (17)  becomes: 


d/i 

d7 


The  energy  moment  of  the  Boltzmann  equation  can  be  written  as  (4:1 1 1) 


(18) 


dn<e>  _  - 

— - +  V  •  n  <  ve  >  -F  •  n  <  v  > 

dt 


<d/n£>coU 


(19) 
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where  F  is  the  force  on  the  electron.  The  second  term  in  this  equation  can  be  zeroed  due  to 
homogeneity.  The  first  term  can  be  expressed  as  two  terms: 


—n <e>  = 

Bt 


B  <  e  >  Bn 

n — - —  +  <£>  — 
Bt  Bt 


(20) 


The  third  term  in  Eq  (19)  can  be  written  as 


Fn*<v>  =  eEn*<\>  =  J*E 


(21) 


Using  Eqs  (18),  (20),  and  (21)  the  er-irgy  moment  of  the  Boltzmann  equation  takes  the  form 

J-E  =  McoW  +  <£>n  (vion  -  vlols)  -  n  (22) 

When  a  self-similar  form  of  the  solution  to  the  Boltzmann  equation  has  been  achieved 
through  the  iteration  method  previously  outlined,  the  last  term  above  will  vanish.  Thus,  the 
power  going  into  the  discharge,  J  •  E,  will  be  balanced  by  the  sum  of  the  inelastic  and  momen¬ 
tum  transfer  losses  contained  in  the  collision  term  (Moo,,),  along  with  the  second  term  on  the  right 
hand  side.  This  term  represents  the  energy  lost  in  bringing  the  secondary  electrons  produced  in 
ionization  up  to  the  average  energy  of  the  distribution,  as  well  as  a  correction  factor  caused  by 
the  loss  of  electrons  through  processes  such  as  attachment,  diffusion,  and  recombination. 

Eq  (22)  proves  to  be  very  useful  in  describing  the  conservation  of  energy  in  the  plasma  dis¬ 
charge.  This  principle  will  serve  as  a  chief  diagnostic  tool  in  determining  whether  each  of  the 
electron  interactions  considered  is  correctly  incorporated  into  the  solution  method. 
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External  Sources 


In  order  to  study  the  plasma  within  a  hydrogen  multicusp  discharge,  excimer  laser,  or  electron 
beam  assisted  discharge,  the  external  source  of  electrons  must  be  included  into  the  solution 
method  of  the  Boltzmann  equation.  Following  closely  the  work  of  Rockwood,  A.E.  Greene  and 
C.J.  Elliot  (6)  developed  a  very  similar  method  of  solution  for  the  specific  application  of  the 
electron  beam  pumped  plasma  of  an  excimer  laser.  In  their  development  of  the  equations  that 
describe  the  physics  involved,  they  included  the  source  term  of  the  e-  beam,  as  well  as  any  loss 
terms  that  were  considered  (in  this  case,  recombination  was  considered  to  be  the  dominate  loss). 
In  order  to  include  these  terms,  reconsider  Eq  (3)  with  a  non-zero  source  and  loss  term.  Under 
such  conditions,  Eq  (23)  is  iteratively  solved  until  convergence. 

n(f+Af)  =  (I-CA/)_1(n(0  +  SA/)  (23) 

The  source  term  here  represents  the  value  of  the  source  at  each  of  the  energy  bins  along  the 
energy  axis.  For  an  electron  beam  with  a  sufficiently  narrow  energy  spread,  this  vector  will  con¬ 
sist  of  zeroes  everywhere  except  at  the  energy  bin  containing  the  e-  beam. 

In  modeling  a  magnetic  multicusp  discharge  in  hydrogen,  Bretagne  (3),  et  al,  has  developed  a 
code  that  also  models  a  source  term  within  a  plasma.  Their  method  of  solution  again  follows  that 
of  Rockwood,  with  the  inclusion  of  the  source  term  and  a  loss  term.  The  dominate  electron  loss 
processes  in  this  case  were  assumed  to  be  recombination  and  diffusion  to  the  walls  of  the  device. 
The  solutions  obtained  by  Bretagne  were  for  the  inclusion  of  an  electron  source  at  90  eV  in 
hydrogen.  This  is  high  when  the  mean  electron  energy  is  of  the  order  of  5-10  eV. 

In  order  to  carry  out  the  calculation  of  the  distribution  of  electrons  to  these  higher  energies, 
Bretagne  adopted  a  non-uniform  differential  energy  element.  In  this  scheme,  the  energy  widths 
were  allocated  logarithmically,  increasing  with  increasing  energy.  Using  this  method,  it  is 
assumed  that  the  distribution  changes  most  rapidly  at  the  lower  energies  where  the  energy  bin 
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widths  arc  small,  and  less  so  at  the  higher  energies  where  the  bin  widths  are  large.  Ideally,  what 
is  required  in  a  solution  method  is  the  allocation  of  a  denser  mesh  of  energy  bins  in  the  vicinity 
of  inelastic  process,  since  this  is  where  the  electron  distribution  will  change  most  rapidly.  This  is 
due  to  the  depleting  effect  that  such  a  process  imposes  on  the  distribution  as  a  result  of  the 
energy  loss  associated  with  the  inelastic  process.  The  energy  loss  of  an  electron  involved  in  such 
a  collision  would  be  the  threshold  energy  of  the  associated  cross  section. 

While  Bretagne’s  method  could  work  well  for  hydrogen,  in  which  most  of  the  inelastic  pro¬ 
cesses  occur  at  energies  which  are  relatively  low  (0.5  -  12  eV)  compared  to  the  energy  of  the 
source,  it  would  seem  to  lack  the  flexibility  needed  to  study  other  gases,  such  as  the  rare  gases,  in 
which  inelastic  processes  can  occur  at  energies  of  20  eV  or  more  (as  in  the  case  of  Helium).  In 
such  a  case,  it  may  be  possible  for  the  energy  bins  to  be  too  wide  at  these  higher  energies, 
thereby  rendering  invalid  the  basic  assumption  used  in  the  finite  differenced  approach. 
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III.  Theory 


Having  reviewed  the  theory  behind  the  basic  Boltzmann  equation  and  the  computational 
method  utilized  in  determining  its  solution,  the  specific  requirements  of  the  plasma  devices  under 
consideration  need  to  be  addressed.  These  items  are  the  processes  of  electron  attachment  and 
electron  -  electron  collisions  along  a  non-uniform  energy  grid,  electron  beams,  and  methods  of 
dynamically  allocating  a  non-uniform  energy  grid.  In  addition,  the  theory  of  steady  state  nega¬ 
tive  mobility  is  addressed. 


Attachment 

The  process  of  electron  attachment  is  very  important  in  the  study  of  certain  plasma  devices, 
particularly  excimer  lasers.  These  devices  are  at  the  forefront  of  technology  in  the  efficient  pro¬ 
duction  of  UV  and  near  UV  power  (18:346).  In  electron  beam  pumped  excimers,  it  is  the  attach¬ 
ment  of  electrons  to  the  halogen  that  forms  the  negative  ion  used  in  the  reaction  process 
(18:348): 


Ar+  +  Fr  +  Ar  (Ar+F')*  + Ar  (24) 

Thus,  in  e-  beam  pumped  excimers,  the  process  of  attachment  serves  as  an  essential  ingredient  in 
the  production  of  the  excited  state.  However,  this  is  not  the  case  in  discharge  pumping  (18:349): 


Ar*  +  F2  -4  ( Ar*F)  +  F 


(25) 


20 


Here,  the  attachment  process  serves,  at  best,  only  as  an  auxiliary  pumping  channel.  The  benefits 
of  the  production  of  ionic  Flourine  are  countered  by  the  detachment  process,  which  can  lead  to 
an  unstable  discharge,  as  pointed  out  by  Verdeyen  (18:349).  So  then,  no  matter  what  the  pump¬ 
ing  mechanism,  electron  attachment  is  an  essential  consideration  in  the  analysis  of  the  excimer 
laser. 

The  attachment  process  is  taken  into  account  by  adding  to  the  collision  term  of  Eq  (4)  the 
attachment  rate: 


^rUch  =  -n<Nv(e)allI14Ch(e)  (26) 

This  attachment  rate  can  be  easily  incorporated  into  the  finite  differenced  method  of  solution  as 
pan  of  the  C  matrix  of  Eq  (9).  The  loss  of  electrons  from  the  system  involves  a  corresponding 
loss  of  energy.  This  energy  loss,  unlike  those  associated  with  ionization  and  excitation,  is  not 
associated  with  a  fixed  threshold  value.  Rather,  the  energy  loss  in  the  attachment  process  can  be 
any  energy  within  the  domain  of  the  attachment  cross-section.  The  electron  that  undergoes 
attachment  in  energy  bin  k  will  lose  an  amount  of  energy  equal  to  the  center  of  bin  energy  for  the 
k*  bin.  Additionally,  in  order  to  achieve  an  energy  balance  once  the  steady-state  (or  self-similar) 
form  of  the  solution  has  been  established,  a  correction  term  for  the  electron  loss  must  be  taken 
into  account,  as  shown  in  Eq  (22). 

Eq  (22)  can  be  rewritten  slightly,  showing  explicitly  the  balance  between  the  energy  into  the 
discharge  and  the  various  inelastic  processes,  along  with  momentum  transfer  losses: 

J-E  =  MT+«ei„,vi„>+<e>v1J  + 

(<eio»Vl<,.>-<E>  VkJ  +  <  ^exc^exc  ^  ^  ^exc^iuper  ^  (27) 
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Investigation  into  Eq  (27)  reveals  several  interesting  cases.  If  the  attachment  frequency  is 
constant,  then  the  average  energy  of  the  electron  distribution  is  exactly  the  same  as  if  there  were 
no  attachment  present.  The  presence  of  the  constant  attachment  collision  frequency  would 
deplete  a  constant  percentage  of  electrons  from  each  energy,  thus  the  form  of  the  distribution 
would  remain  the  same  as  if  attachment  were  not  present.  The  total  number  of  electrons  would, 
of  course,  be  less  with  an  attachment  process. 

Consider  the  case  of  an  attachment  frequency  which  increases  with  energy.  Here,  the  attach¬ 
ing  process  would  favor  the  high  energy  electrons  over  the  slower  ones,  thus  depleting  the  disui- 
bution  at  the  higher  energies.  The  result  would  be  an  EEDF  with  a  lower  average  energy. 

Finally,  consider  an  attaching  frequency  which  decreases  with  energy.  In  this  case,  the 
attachment  process  would  tend  to  deplete  electrons  at  the  lower  energies,  resulting  in  a  distribu¬ 
tion  which  would  have  a  higher  mean  energy  than  otherwise. 

Steady  State  Negative  Mobility 

It  is  this  last  attachment  case  which  forms  the  origin  for  the  possibility  of  a  plasma  discharge 
with  a  negative  steady  state  total  electron  mobility.  In  general,  the  mobility  can  be  expressed  as: 

V, 

n  =  y  (28) 

where  E  is  the  externally  applied  electric  field,  and  vd  is  the  drift  velocity.  Using  the  classical 
two  term  spherical  harmonic  expansion  with  Eq  (13)  for  the  drift  velocity,  along  with  /,  given  by 
(10:471): 
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(29) 


Ed/,  1 
N  de  o„  +  L  cA 


A 

where  oh  is  the  inelastic  cross  section  for  process  h,  the  mobility  can  be  rewritten  as  (15:1594): 


=  -fVTq/ 


d/0(e) 


Na„(e)  d£ 


de 


(30) 


It  has  been  assumed  that  the  sum  of  the  inelastic  cross  sections  is  much  smaller  than  the  cross 
section  for  momentum  transfer.  Eq  (30)  can  be  integrated  by  parts,  with  the  integral  term 
becoming: 


fe3/„  _  tf0  -  r  [  1 

J  am  3e  cm  o 


efo 


3(0 


oi  3e 


kie 


(31) 


J 


0  0 
The  first  term  on  the  right  hand  side  is  zero.  The  evaluation  of  this  term  at  infinity  is  zero  due  to 
the  nature  of  the  distribution  function  in  that  limit,  and  the  evaluation  at  the  lower  limit  is  easily 
handled  remembering  that  the  cross  section  is  finite  at  zero  energy.  Thus,  a  criteria  for  negative 
mobility  is  established: 


dam  a(e) 
de  e 


(32) 


Graphically,  this  criteria  can  be  envisioned  by  drawing  a  line  from  the  origin  to  the  g„(e)  at  a 


particular  energy  (Figure  1).  If  the  slope  of  the  momentum  transfer  cross  section  at  this  energy  is 
greater  than  the  slope  of  the  line  drawn,  then  the  criteria  for  negative  mobility  is  met.  It  is  seen 
that  this  criteria  is  met  by  Argon  in  the  energy  range  of  approximately  0.3  to  10.0  eV. 
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Figure  1.  At  momentum  transfer  cross  section  (21 : 1543) 


An  excellent  description  of  the  physical  process  at  work  in  a  discharge  with  a  steady  state 
negative  mobility  is  given  by  Rozenburg,  et  al  (15:1594).  Briefly,  however,  it  is  the  presence  of 
an  attaching  cross  section  at  low  energies  which  can  lead  to  the  effect.  I?  the  electron  distribu¬ 
tion  is  shaped  such  that  the  majority  of  the  integration  of  Eq  (30)  is  carried  out  in  the  energy 
range  over  which  Eq  (32)  is  valid,  then  the  mobility  will  be  negative.  Electrons  travelling 
against  the  direction  of  the  external  electric  field  will  tend  to  gain  energy  from  the  field.  How¬ 
ever,  this  energy  gain  will  be  transfered  to  the  gas  molecules  through  elastic  collisions  rather 
quickly,  because  the  momentum  transfer  cross  section  increases  rapidly  with  energy.  Thus,  these 
electrons  will  eventually  end  up  with  a  velocity  component  in  the  direction  of  the  field,  at  low 
fields.  The  electrons  originally  travelling  in  the  same  direction  as  the  electric  field  will  slow 
down  due  to  the  force  of  the  field  and  will  lose  energy.  If  a  large  cross  section  for  attachment  is 
present  at  the  lower  energies,  then  there  will  be  a  large  probability  that  the  electron  will  be  lost. 


24 


Additionally,  since  the  electrons  are  slowing  down,  they  will  undergo  fewer  and  fewer  elastic 
collisions  with  the  gas  molecules,  thus  they  will  not  tend  to  assume  a  spherical  velocity  distribu¬ 
tion.  These  two  processes  will  compete  with  one  another,  and  conditions  may  be  reached 
wherein  the  electrons  will  have  a  mean  velocity  component  in  the  direction  of  the  electric  field, 
and  therefore  have  a  negative  mobility. 

Electron-Electron  Interactions 

The  effect  of  electron-electron  collisions  in  a  plasma  is  to  drive  the  distribution  toward  a 
Maxwellian  (4:66).  This  effect  is  particularly  important  for  discharges  with  high  fractional  ion¬ 
ization,  and  low  E/N  values  (14:2350).  Unlike  the  inelastic  collisions  that  are  incorporated  into 
the  C  matrix  of  Eq  (9),  the  electron-electron  interactions  are  nonlinear  and  therefore 
time-dependent.  Being  described  by  the  coulomb  collision  integral,  the  effect  of  a  collision 
between  electrons  depends  on  the  entire  distribution  of  electrons  in  energy  space.  This  nonlinear 
interaction  is  incorporated  into  the  solution  method  by  casting  this  collision  integral  into  the 
energy  space  flux  divergent  formalism.  In  doing  so,  the  two  physical  properties  of  conservation 
of  particles  and  energy  must  be  ensured.  Additionally,  the  electron  -  electron  collisions  must 
tend  to  drive  the  distribution  to  a  Maxwellian  over  time,  reaching  a  Maxwellian  in  the  steady 
state  when  only  electron  interactions  are  considered. 

The  recasting  of  the  collision  integral  into  the  flux  divergent  formalism  is  accomplished  by 
the  following  steps,  and  is  developed  in  detail  in  Appendix  B.  Eq  (7),  which  represents  the 
Fokker-Planck  approximation  to  the  elec  iron -electron  collision  integral,  is  finite  differenced 
using  a  simple  backward  difference  scheme.  In  doing  so,  the  effect  of  the  non-uniform  energy 
axis  is  taken  into  account  by  defining  the  average  electron  number  density  between  two  adjacent 
bins  as 
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AM*  +  i  +  A£*  +  1nA 


(33) 


nk  +  \n 


Ae*  +  A £*+1 


while  the  term  representing  the  energy  space  divergence  of  the  electron  number  density  becomes 


dfljfc  + 1/2  _  nk  +  \~nk 

£*  + 1  “ 


(34) 


Notice  that  with  constant  energy  bin  widths,  these  representations  will  recover  the  usual  expres¬ 
sions  for  the  given  quantities. 

Next,  the  integrals  contained  in  Eq  (7)  are  replaced  by  summations  over  the  entire  energy 
range,  and  like  terms  of  the  electron  number  density  are  grouped  together  and  expressed  in  the 
same  formalism  used  previously  to  describe  the  field  and  elastic  driven  fluxes  found  in  Eq  (8). 
Eq  (35)  and  Eq  (36)  result  Bt  +  1 ,,  the  coefficient  of  ni  +  ),  represents  the  downflux  of  electrons 
going  from  energy  bin  k+1  to  k  by  electrons  going  from  bin  1  to  1+1.  The  coefficient  of  n*_,, 
which  is  A*_,  „  represents  the  rate  of  upflux  of  electrons  from  k-1  to  k,  by  electrons  going  from  1 
to  1-1. 


dnk 

~dT=a 


+  b  k  +  \nk  +  \  (a  k  +  b  k)nk 


where  the  coefficients  are  expressed  as: 


(35) 


a  *-i  ~  £ \-uni  b  k+i  -  S  (36) 

In  order  to  ensure  particle  conservation,  boundary  conditions  are  applied  to  Eq  (35)  such  that 
particles  cannot  be  lost  from  the  energy  grid  during  the  iteration  process  used  to  determine  the 
solution.  This  condition,  which  requires  that  electrons  cannot  be  demoted  from  the  lowest 
energy  bin,  nor  can  they  be  excited  from  the  highest  energy  bin,  is  expressed  mathematically  as: 
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An.;  =  Bu  =  0  for  all  /  £  [/,A/] 


(37) 


Energy  conservation  requires  that 


oo 

f  d  N 

£n(e)d£  =  —  X  tnk  Ae.  =  0 
J  dt*  =  i 


Carrying  out  the  algebraic  operations  required  in  Eq  (38)  results  in 


f  Ae;  +  1  +  A£,  ] 

B‘-'  =  (39: 

In  the  steady  state  with  a  Maxwellian  distribution,  the  promotion  and  demotion  interpretation 
for  the  Ak  ,’s  and  B*  ,’s  leads  to  the  following  balance  between  the  two: 
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This  expression  states  that  the  rate  of  electron  promotion  from  energy  bin  k  to  k+1  by  elec¬ 
trons  going  from  bin  /  to  /  -  1  is  balanced  exactly  by  the  rate  of  demotion  of  electrons  from  bin 
k+1  to  k  by  electrons  going  from  /  -  1  to  /.  Using  Eq  (38)  and  the  expression  for  a  Maxwellian 
energy  distribution,  Eq  (40)  can  be  expressed  as: 
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where  T  is  the  electron  temperature.  The  proper  construction  of  the  A  matrix  will  ensure  that  the 
distribution  will  be  a  Maxwellian  in  the  steady  state.  The  use  of  Eq  (40)  guarantees  that  energy 
will  be  conserved  in  the  closed  system,  while  the  boundary  conditions  defined  by  Eq  (38)  will 
ensure  the  conservation  of  electrons. 


Variable  Energy  Grid 


The  variable  energy  grid  is  envisioned  as  not  only  a  method  to  increase  the  energy  range  of  a 
calculation,  but  also  as  a  means  of  increasing  the  accuracy  of  transport  parameters  and  kinetic 
rates  over  the  uniform  grid.  By  increasing  the  number  of  energy  bins,  while  holding  the  maxi¬ 
mum  energy  constant,  the  finite  differenced  representation  of  the  flux  terms  becomes  a  better 
approximation  to  reality,  thus  increasing  the  accuracy  of  the  computed  distribution  to  the  actual 
distribution.  In  so  doing,  it  is  expected  that  the  resultant  transport  parameters  and  kinetic  rates 
computed  from  the  distribution  would  also  increase  in  accuracy.  This  trend  has  a  practical  limit 
however.  As  the  width  of  the  energy  bin  becomes  sufficiently  small,  the  computer  roundoff 
error  will  begin  to  accumulate  and  result  in  significant  deviations  from  the  true  values. 

When  considering  a  non-uniform  energy  grid,  the  question  that  naturally  arises  is  how  this 
mesh  should  be  allocated.  In  the  electron  beam  sustained  plasma  discharges,  the  logarithmically 
allocated  axis  has  been  used  by  Bretagne,  which  seems  to  be  a  useful  method  for  H2  and  perhaps 
other  molecular  gases.  However,  it  may  not  be  the  logical  choice  for  gases  with  inelastic  pro¬ 
cesses  occurring  at  higher  energies  than  those  present  in  Hz  such  as  the  rare  gases.  This  is  espe¬ 
cially  true  if  the  ionization  and  excitation  rates  of  such  gases  are  of  interest.  Additionally,  this 
method  will  indiscriminately  allocate  an  energy  axis  with  no  regard  as  to  the  inelastic  cross 
sections  present.  It  would  seem  that  the  logarithmic  method  is  used  primarily  for  its  ease  of 
implementation. 
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A  better  method  of  allocation  would  logically  take  the  inelastic  cross  sections  and  threshold 
energies  into  account.  These  inelastic  processes  would  decrease  the  distribution  very  effectively 
at  the  threshold  energies  and  beyond,  thus  the  increased  bin  density  at  these  energies  would  serve 
to  increase  the  accuracy  of  the  finite  differenced  representation  of  the  differential  equation. 

Since  average  energy  and  drift  velocity  are  parameters  that  are  mostly  dependent  on  the  bulk  of 
the  distribution,  it  is  essential  that  the  mesh  used  accurately  calculates  this  region  of  the  EDF. 
Additionally,  excitation  and  ionization  rates  arc  parameters  that  mostly  depend  on  the  tail  of  the 
distribution,  therefore  it  is  important  to  be  able  to  accurately  calculate  this  portion  of  the  EDF  as 
well.  The  best  possible  mesh  to  calculate  the  former  is  clearly  not  the  best  to  calculate  the  latter. 
The  challenge  is  to  be  able  to  calculate  both  regions  of  the  distribution  with  greater  accuracy  than 
a  uniform  mesh  of  equal  grid  points. 

One  method  that  can  be  used  to  define  an  energy  grid  has  been  used  by  Nickel  (1 1:18)  to 
calculate  the  radiation  from  various  spectral  lines.  In  his  method,  a  function  F  is  defined  by  the 
relation 


J  G(e)de 

He)  =  ~ -  (42) 

J  G(e)de 

0 

where  G  is  some  function  yet  to  be  determined,  but  which  depends  in  this  instance  on  the  cross 
sections  of  the  elastic  and  inelastic  processes.  Notice  that  the  value  of  F  can  range  between  the 
limits  of  0  and  1.  The  energy  axis  is  divided  into  N  equally  spaced  points,  and  the  value  of  the 
function  F  is  computed  at  each  point  These  points  are  then  graphed  versus  energy  as  shown  in 
Figure  2,  where  a  10  point  uniform  grid  has  been  mapped  onto  a  10  point  variable  grid  by  using 
the  Maxwellian  distribution  as  the  function  G.  The  y  axis,  which  ranges  from  0  to  1,  is  also 
divided  into  N  equally  spaced  points.  Each  of  these  points  is  then  projected  horizontally  onto  the 
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graph  of  F  vs  energy  until  intersection,  and  then  down  vertically  to  the  energy  axis.  These  inter¬ 
sections  along  the  energy  axis  now  define  the  new  energy  axis  along  which  the  solution  to  the 
Boltzmann  equation  will  be  computed.  It  is  clearly  ieen  from  Eq  (42)  and  figure  2  that  in  energy 
regions  where  the  function  G,  and  consequently  function  F,  changes  rapidly,  the  mesh  of  energy 
bins  will  increase.  On  the  other  hand,  in  regions  where  F  changes  little,  few  energy  bins  will  be 
allocated.  This  is  precisely  the  allocation  scheme  that  is  required,  the  only  question  being  what 
does  function  G  look  like?  Notice  that  if  the  function  G  were  defined  such  that  F  was  a  straight 
line,  then  the  mapping  would  be  uniform  to  uniform. 


Figure  2.  Maxwellian  energy  distribution  with  electron  temperature  of  2  units  and  the  corre¬ 
sponding  function  F,  where  G  is  given  by  the  number  distribution 

The  choices  of  the  function  G  can  range  from  the  electron  distribution  itself  to  some  moment 
of  the  distribution,  such  as  the  energy  moment  Unfortunately,  these  energy  bin  allocation  func¬ 
tions  require  the  distribution  function  itself,  and  this  is  the  solution  of  the  Boltzmann  equation. 
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Thus  these  allocation  choices  require  either  first  solving  the  Boltzmann  equation  on  a  uniform 
energy  grid,  and  then  re-solving  it  on  the  non-uniform  grid,  or  accepting  some  assumptions  that 
would  allow  the  calculation  of  an  approximation  to  the  electron  distribution  quickly,  without  the 
need  to  iterate  until  convergence.  The  former  choice  is  not  ideal,  due  to  the  time  requirements. 
Thus  the  latter  alternative  becomes  the  more  attractive  method. 

Using  the  commonly  accepted  two  term  spherical  harmonic  expansion,  the  electron  distribu¬ 
tion  can  be  represented  a sf=f0  +/, cos(9),  where  f0  and  /,  are  the  solutions  of  a  pair  of  coupled 
differential  equations  (10:471): 


E/N  9 
3e  de 


'V'' 


N  de  + 


(a+io,>,= 


0 


(43a) 

(43b) 


where  Q  is  a  cross  section  and  e'  =  e+ 8e*.  The  subscript  h  denotes  an  inelastic  process,  while 

the  subscript  r  denotes  weighting  each  momentum  transfer  cross  section  by  the  mole  fraction  of 
the  gas.  The  subscript  m  denotes  the  momentum  transfer  cross  section  weighted  by  mole  frac¬ 
tion  and  molecular  weight  for  the  gas.  Following  Holstein’s  example  (8:372),  Eq  (43a)  is  inte¬ 
grated  with  respect  to  energy,  then  Eq  (43b)  is  substituted  into  Eq  (43a),  and  again  integrated 
with  respect  to  energy.  The  pair  of  coupled  differential  equations  now  becomes  a  single  implicit 
integral  equation: 
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(44) 


/„(e)=eJtp[-^;(^  j*Q,WQ,00d* 

f  Q,W  '"f* 

+J^w  j  Q‘(y)/-(y)ydydxI 

where  Q,  is  the  sum  of  all  cross  sections.  Now,  if  some  simplifying  assumption  can  be  made  to 

eliminate  f0  on  the  right  hand  side  of  Eq  (44),  then  this  function  could  be  calculated  quickly  and 
may  serve  as  a  good  choice  to  be  used  in  Eq  (42)  in  defining  a  non-uniform  energy  mesh.  If  it  is 
assumed  that  fQ  remains  relatively  constant  over  the  interval  £  to  E  + Ae*,  then  this  simplification 
will  occur.  The  result  is  an  equation  for  the  electron  energy  distribution  function  which  is  a 
function  only  of  energy,  E/N,  gas  type,  and  the  cross  sections  for  the  gas.  Since  these  are  all 
parameters  that  are  known  prior  to  the  start  of  the  numerical  calculation  of  the  solution  to  the 
Boltzmann  equation,  this  function  can  be  used  quickly  in  defining  a  new  energy  grid  through  Eq 
(42). 

Once  the  new  energy  grid  has  been  established,  and  the  solution  to  the  Boltzmann  equation 
computed  along  it,  some  criteria  must  be  used  by  which  the  utility  of  the  non-uniform  method 
may  be  evaluated.  Analytic  solutions  to  the  Boltzmann  equation  exist  for  the  special  cases  of 
constant  momentum  transfer  cross  section  and  for  constant  momentum  transfer  collisional  fre¬ 
quency.  In  both  of  these  cases,  any  inelastic  cross  section,  Q,,  has  been  zeroed  out.  Therefore 
the  double  integral  on  the  right  hand  side  of  Eq  (44)  vanishes.  The  remaining  integral  can  be 
handled  analytically,  allowing  calculation  of  such  parameters  as  average  energy  and  drift  veloc¬ 
ity  through  the  use  Eq  (15)  and  Eq  (13)  respectively.  Therefore,  in  these  cases,  the  average 
energy  and  drift  velocity  calculated  using  uniform  and  non-uniform  energy  grids  can  be 
compared  to  the  analytically  determined  ones.  However,  since  the  primary  utility  of  the  non- 
uniform  energy  grid  most  likely  lies  in  the  case  of  when  inelastic  cross  sections  are  included,  the 
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analytic  calculation  of  these  parameters  is  precluded.  In  this  case,  a  good  figure  of  merit  is  a 
calculation  of  the  parameters  carried  out  using  a  uniform  energy  grid  consisting  of  a  sufficient 
number  of  points  such  that  the  calculated  parameters  do  not  change  much. 

Electron  Beams 

Electron  beams  can  be  included  into  the  method  of  solution  by  returning  to  Eq  (9).  If  the 
source  and  loss  terms  in  Eq  (3)  are  non-zero,  then  Eq  (9)  becomes: 

d nk 

— i=XCwn,  +  S-L  (45) 

at  i 

where  the  units  of  the  source  term  S  and  loss  term  L  would  be  cm"V\  Applying  Euler’s  formula 
for  the  derivative  and  solving  for  n(t  +  A/),  the  following  is  obtained: 

n(t  +  At)  =  (I  -  CAf)-1  {(S  -  L)A/  +  n(r )}  (46) 

The  solution  to  the  Boltzmann  equation  is  found  at  each  time  interval  t  +  Ar,  and  the  self- 

similar  solution  is  found  by  iterating  Eq  (46)  until  the  programmed  convergence  criteria  is  met. 

In  analyzing  the  effect  of  only  the  e-  beam  on  the  distribution,  it  is  desirable  to  zero  the  elec¬ 
tric  field.  Doing  so  with  the  finite  differenced  equations  used  in  previous  codes  created  numer¬ 
ical  instabilities  arising  from  the  recoil  flux  term,  J ^  (6:2948).  In  those  representations  for  the 
flux,  the  reduction  of  the  electric  field  below  some  minimum  value  allowed  promotion  and 
demotion  rate  coefficients  (ak’s  and  bk’s)  that  were  almost  all  negative  at  some  values  of  energy. 
With  no  forces  acting  on  the  closed  system  of  electrons,  what  is  required  of  the  electron  energy 
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distribution  is  to  obtain  a  Maxwellian  whose  average  energy  is  l.S  times  the  gas  temperature. 
The  finite  differenced  representation  for  the  recoil  flux  divergent  term  as  achieved  in  Appendix 
A  preserves  this  property,  thus  maintaining  numerical  stability  over  a  range  of  temperatures. 
The  source  term  S  can  be  written  as 


s=£  (47) 

where  V  is  the  plasma  volume,  and  I  is  the  current  at  the  source  energy.  Eq  (47)  is  written  for  a 
point  source  in  energy.  If  a  source  is  used  that  is  broad  in  relation  to  the  energy  bins  around  it, 
then  this  term  would  need  to  be  spread  out  into  several  bins. 

The  addition  of  the  electron  beam  source  also  requires  its  inclusion  into  the  energy  balance 
equation  Eq  (22),  where  now  the  J  •  E  term  is  replaced  (or  augmented)  by  the  power  into  the  dis¬ 
charge  by  the  e-  beam,  namely: 


bum  =  («) 

where  the  units  of  this  power  term  are  eVs_,cm'\  a  power  density. 

The  effect  of  having  a  source  of  electrons  at  a  high  energy  is  to  produce  ionization  and  excita¬ 
tion  of  the  gas  molecules.  These  processes  all  constitute  losses  of  energy  for  the  source  elec¬ 
trons,  thus  they  will  be  displaced  down  to  the  lower  energies  by  an  amount  equal  to  the  threshold 
energy  of  the  process.  If  the  source  is  at  large  enough  energies,  these  electrons  can  be  displaced 
along  the  energy  grid  many  times  before  elastic  collisions  with  the  gas  become  dominate.  The 
interactions  between  the  ionizing  electrons  and  the  secondary  electrons  produced  during  ioniza¬ 
tion  gives  rise  to  what  Bretagne  (1:1119)  calls  energy  sharing.  In  the  present  Boltzmann  code, 
the  secondary  electrons  are  all  assumed  to  appear  at  the  lowest  energy  bin,  while  in  reality  some 
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energy  is  shared  between  it  and  the  ionizing  electron,  resulting  in  some  secondary  electrons 
appearing  at  higher  energies.  This  energy  sharing  tends  to  smooth  the  distribution  when  com¬ 
pared  to  the  classical  treatment  of  ionization  that  neglects  it.  However,  according  to  Bretagne, 
this  classical  treatment  of  the  ionization  process  can  be  used  to  accurately  model  the  electron 
beam  sustained  discharge  as  long  as  1)  low  energy  inelastic  process  are  not  dominate  and  2)  the 
energy  of  the  source  is  not  too  high  ( Ep  <  about  100  eV)  (1:1201).  The  first  condition  can 
readily  be  met  in  the  rare  gases,  however  in  molecular  gases  care  must  be  taken  to  ensure  that 
this  is  true. 


IV.  Computational  Method 


Program  Structure 

The  Boltzmann  code  developed,  Megaboltz,  solves  the  time-dependent  Boltzmann  equation. 
Output  consists  of  the  EEDF,  the  kinetic  rates  for  the  inelastic  processes  considered,  and  the 
transport  parameters.  An  energy  balance,  which  is  the  primary  compututational  diagnostic,  is 
also  a  monitored  output.  A  user’s  guide  to  the  code  is  included  in  Appendix  C,  however  the 
principal  routines  will  be  briefly  outlined  in  this  section. 

The  user  can  select  from  a  menu  regarding  the  type  and  amount  of  each  gas  present  in  the 
plasma,  the  electric  field  to  neutral  number  density  ratio  (E/N),  the  gas  temperature,  pressure, 
and  electron  number  density.  Additionally,  the  user  can  choose  to  see  or  plot  the  interpolated 
cross  sections  for  any  of  the  kinetic  processes,  as  well  as  turn  on/off  the  electron-electron  colli¬ 
sions  and  electron  beam.  The  source  energy  of  the  electron  beam  is  also  set  by  the  user  on  the 
menu.  Finally,  the  energy  grid  is  selected,  whereby  the  number  of  energy  bins  along  the  energy 
axis  and  the  maximum  energy  are  defined.  This  grid  can  be  constructed  in  a  number  of  ways: 
uniformally,  logarithmically,  or  dynamically.  When  the  energy  grid  is  constructed  dynamically, 
the  elastic  and  inelastic  cross  sections  will  have  some  influence  on  the  location  and  density  of  the 
energy  bins. 

The  value  of  the  independent  variable  used  in  the  numerical  integration  of  the  Boltzmann 
equation  is  also  set  in  "input.com".  This  number  represents  the  time  step  of  integration  used  in 
Eq  (1 1)  and  Eq  (12),  and  would  typically  be  at  least  as  small  as  the  time  scale  for  the  dominate 
electron  collisional  process. 
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In  order  to  properly  invoke  the  code,  the  user  must  define  the  physical  situation  of  interest. 
This  is  done  by  editing  the  file  called  "input.com".  This  file  will  be  read  when  the  program  is 
run,  thus  it  must  be  edited  before  starting  the  program.  This  file  defines  the  basic  problem  to  be 
solved. 

After  the  input  file  "input.com"  is  read,  the  code  will  look  for  the  cross  section  files  of  the 
gases  that  were  selected.  These  cross  section  files  are  external  to  the  code,  and  they  must  be 
formatted  according  to  Appendix  D.  These  files  contain  cross  section  vs  energy  data  stored  as 
tables,  along  with  some  other  required  information.  For  each  type  of  inelastic  process,  a  number 
called  NSTATE  is  read,  which  tells  the  code  the  ratio  of  the  gas  molecules  that  are  in  the  lower 
state  of  the  process  the  cross  section  represents  to  the  total  number  of  gas  molecules.  Excitation 
or  ionization  from  the  ground  state  would  have  NSTATE  values  of  1.0.  An  excitation  from  some 
excited  state  to  a  higher  state  would  have  an  NSTATE  value  of  between  0.0  and  1 .0.  By  assign¬ 
ing  each  process  a  separate  NSTATE  value,  any  of  the  inelastic  processes  may  be  turned  on  or 
off.  For  excitation  processes,  an  additional  variable  is  read  which  defines  the  ratio  of  the 
molecules  in  the  upper  excited  state  to  the  number  in  the  lower  state,  which  may  or  may  not  be 
the  ground  state.  This  value,  NSTAR,  is  used  to  compute  the  superelastic  rates. 

Once  all  of  the  gas  cross  sections  are  read  in,  the  energy  axis  is  defined  according  to  the 
user’s  choice  in  "inpuLcom",  and  the  elastic  and  inelastic  cross  sections  are  interpolated  from  the 
data  tables  at  each  energy  point.  Next,  the  matrix  C  as  used  in  Eq  (9)  is  constructed.  The  ele¬ 
ments  of  this  matrix  include  all  excitation,  ionization,  superelastic  and  attachment  processes. 
Each  element  has  units  of  s'1,  and  is  determined  by  how  the  particular  process  effects  the  flow  of 
electrons  in  energy  space  at  each  bin  energy.  Since  the  number  density  of  electrons,  nk,  is  a 
quantity  that  is  defined  at  the  center  of  energy  bin  k,  the  coefficients  of  nk  are  also  constructed  to 
be  center  of  bin  values. 
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Impact  excitation  and  ionization  will  cause  an  electron  with  energy  £  >  to  lose  an 

amount  of  energy  equal  to  this  threshold  value.  Thus,  an  electron  will  be  lost  out  of  the  higher 
energy  bin,  reappearing  at  the  bin  corresponding  to  e  -  Ae^^,  as  outlined  by  Eq  (49). 


Qu+nij  —  2rNfNjO^£k+m.jv^£k+in.j 

=  -NINJo(i1,)V(ik)  (49) 

where  s  is  the  gas  type,  j  is  the  process,  k  represents  the  energy  bin,  nij  is  the  bin  offset  corre¬ 
sponding  to  a  threshold  energy  of  A£j.  Nj  is  the  NSTATE  value  for  process  j,  and  N,  is  the  num¬ 
ber  density  of  gas  molecules  for  gas  type  s.  The  center  of  bin  energy  value  for  bin  k  is  expressed 
as  £*.  The  first  term  on  the  right  hand  side  represents  gains  into  bin  k  due  to  electrons 
undergoing  collisions  at  energy  k+m,  while  the  second  term  represents  excitation  losses  out  of 
bin  k. 

Superelastic  collisions  are  handled  similarly,  except  now  the  addition  of  electrons  at  energy  e 

are  from  energy  e  -  A£j.  The  electron  undergoing  a  superelastic  collision  at  energy  e  will  be  sent 
to  £  + A£j.  The  superelastic  matrix  elements  are  assigned  by 


=  NN, 


{  £k 


Ek-AEj 


ktgvd^-Aep 


J 


=  -N*N  I 


b^+Af-pVCtg 


(50) 


V 


J 


where  k  -  mj  is  the  energy  bin  corresponding  to  -  AEj,  and  detailed  balance  has  been  used  to 


express  the  superelastic  cross  section  in  terms  of  the  excitation  cross  section. 
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In  the  process  of  attachment,  the  electron  is  lost  from  the  system  -  they  do  not  reappear  any¬ 
where  so  far  as  the  energy  distribution  is  concerned.  This  loss  of  electrons  is  assigned  to  ele¬ 
ments  in  C  by 


Cu,  =  -N,N,utta(ek)V(ek)  (51) 

Once  all  of  the  processes  have  been  loaded  into  the  C  matrix,  the  inverse  of  the  (I-Cdt)  is 
calculated.  In  the  present  version  of  the  code,  this  matrix  inversion  is  carried  out  by  the  IMSL  v. 
9.2  library  routine  LUDATF.  This  routine  performs  the  matrix  inversion  usinfe  ihe  L-U  decom¬ 
position  method  (12:31).  As  is  reported  by  Seger  (16:16)  and  Honey  (9:1 18),  L-U  decomposi¬ 
tion  is  generally  preferable  to  other  methods  such  as  Gauss-Jordon,  Gauss-Seidel,  and  successive 
over  relaxation  (SOR)  in  terms  of  both  speed  and  accuracy.  In  the  L-U  decomposition  method, 
the  inverted  matrix  is  actually  two  matrices,  one  upper  triangular  (U)  and  the  other  lower  triangu¬ 
lar  (L),  stored  as  a  single  matrix.  The  new  /T(f  +  Ar)  vector  must  be  computed  by  forward  substi¬ 
tution  with  the  L  matrix,  and  then  back-solving  with  the  U  matrix.  Once  the  (I  -  CAr)  matrix  is 
triangularized,  which  only  needs  to  be  done  once,  then  the  forward/backward  substitution  is 
carried  out  fairly  quickly. 

Each  time  the  new  n(t  -  A/)  vector  is  computed,  it  is  tested  for  convergence  with  the  previous 

vector,  n(t).  This  convergence  test  is  accomplished  at  each  energy  value,  and  each  number  den¬ 
sity  vector  is  normalized  to  its  respective  total  electron  number  density.  Once  the  convergence 
test  is  passed,  or  if  the  maximum  number  of  iterations  is  reached,  the  last  n{t  +  Ar)  vector 
computed  is  passed  to  a  routine  which  calculates  the  excitation,  ionization,  superelastic,  and 
attachment  rates.  These  rates  are  calculated  by 
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where  the  units  on  the  rates  are  cm3/s. 


(52) 


Energy  Balance 


When  the  rates  have  been  calculated,  the  final  electron  distribution  is  passed  off  to  the  energy 
balance  routine.  This  routine  computes  the  energy  into  the  discharge,  and  the  energy  losses  from 
all  elastic  and  inelastic  processes.  The  balance  of  energy  inputs  with  loss  mechanisms  serves  as 
a  primary  computational  diagnostic  for  the  acceptance  of  any  data  run.  The  energy  sources  con¬ 
sist  of  the  electric  field  and/or  the  electron  beam.  The  electric  field’s  energy  contribution  into 
the  discharge  is  just  J  •  E.  This  can  be  written  as 

^  (53) 

where  ak  represents  the  flux  of  electrons  through  energy  space  from  energy  e*  to  energy  e*+,  due 

to  the  applied  field,  and  defined  as  in  Appendix  A.  Similarly,  bt  represents  the  flux  of  electrons 
from  e*  to  e*  _ ,  due  to  the  applied  field.  The  separate  Ae*’s  illustrate  the  effect  of  the  non- 
uniform  energy  axis.  The  first  term  reflects  the  energy  gained  by  electrons  which  have  velocity 
components  anti-parallel  to  the  electric  field,  therefore  accelerating  to  higher  velocities.  The  bk 
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term  reflects  the  slowing  down  of  the  electrons  which  have  some  velocity  component  in  the  same 
direction  as  the  field.  If  more  electrons  speed  up  than  slow  down,  then  the  ,>'d  nave  a  net 
positive  input  of  power  into  the  discharge.  If  more  electrons  slow  down,  then  the  field  will  gain 
energy  from  the  electrons.  This  last  case  is  predicted  from  the  effects  of  a  attachment  cross 
section  at  low  energy,  and  a  momentum  transfer  cross  section  that  increases  with  energy 
(15:1596).  The  energy  into  the  discharge  is  given  by  Eq  (53). 

The  energy  losses  consist  of  elastic  and  inelastic  collisions.  The  rate  of  energy  lost  in  the 
elastic  collisions  can  be  written  as 

lei  =  -Z[(ak-ak)A£.k-(bk-bk)A£k_x]nk  (54) 

The  excitation  loss  is  expressed  as 

^|exc  =  -lk  KsNjnkOj(lk)V(ik)A£SJ  (55) 

where  Ae  is  the  excitation  energy  threshold  for  process  j  of  gas  s  and  N,  is  the  ratio  of  the  number 
density  of  gas  molecules  in  the  lower  state  of  process  j  to  the  number  density  in  the  ground  state. 
Thus,  for  ground  state  excitation  this  number  would  be  1 .0,  whereas  for  excitation  from  some 
excited  state  it  would  less  that  1.0.  The  ionization  energy  loss  is  similar  to  the  excitation  loss, 
with  the  additional  term  reflecting  the  energy  lost  in  bringing  all  the  secondary  electrons  up  to 
the  average  energy  of  the  distribution.  It  is  expressed  as 
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-J-L  =  -JtN,N,«(aj(e,)V(et)A£,J 

Xe4n4 

-V-  I  N,N;a/I4)V(i4)n4  (56) 

S,j,k 

k 

where  the  second  rate  sum  is  carried  out  over  all  ionization  processes  j. 

The  superelastic  loss  is  really  an  energy  gain,  and  can  be  written  as 


~r  L,cl,  =  -  I  N,N,«to(e,)V(E4)£4 

at  s,j,k 

Xe4n4 

+V-  I  N,NAa/i4)V(i4)  (58) 

Z*nk  s,j,k 
k 

where  the  second  rate  sum  is  carried  out  over  all  attachment  processes  j.  The  second  term  in  Eq 
(58)  represents  the  correction  to  the  energy  lost  in  attachment  due  to  loss  of  electrons,  as  Eq  (22) 
suggests. 
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Transport  Parameters 


The  calculation  of  the  transport  parameters  is  accomplished  after  MEGABOLTZ  computes 
the  energy  balance.  These  parameters  are  drift  velocity,  average  energy,  characteristic  energy, 
and  the  free  diffusion  coefficient.  The  drift  velocity  Eq  (13)  can  be  expressed  in  terms  of  the  rate 
of  energy  gained  by  the  electrons  from  the  electric  field  (14:2350): 


vd  = 


E  In* 

k 


(59) 


where  Eg  was  calculated  in  Eq  (53)  as  part  of  the  energy  balance.  The  diffusion  coefficient  is 


calculated  with  Eq  (14)  as 


D'  ■  5 


\™J  *  IqsOs(Ek) 


where  fk  =  - 


ru. 


e*  Ae*  X  nk 

k 


(60) 


Once  the  drift  velocity  and  diffusion  coefficient  are  calculated,  then  the  characteristic  energy 
can  be  calculated  by  Eq  (16),  which  can  be  expressed  as 


ED, 


E*  = 


(61) 


since  p  =  vd/E  . 

The  average  energy  is  calculated  by 
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<  E  > 


(62) 


Im* 

k _ 

I  nk 

k 

Once  these  parameters  arc  calculated,  they  can  be  compared  to  experimental  data.  Very 
often,  it  is  through  just  such  a  comparison  that  a  particular  gas’  cross  section  that  s  unknown  in 
some  energy  range  is  refined  or  even  defined  altogether.  It  is  these  transport  parameters  that  link 
the  distribution  function  to  observables  in  the  real  world. 
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Validation  with  Non-Uniform  Energy  Grid 


The  extension  of  the  energy  range  of  the  solution  is  critical  in  being  able  to  include  sources  of 
electrons  in  the  computation  of  the  distribution  function.  This  energy  extension  involved  re¬ 
deriving  the  field  and  recoil  electron  flux  terms  for  a  non-uniform  energy  grid.  Several  methods 
were  used  in  order  to  validate  the  Boltzmann  code  with  the  non-uniform  energy  mesh.  First, 
calculations  of  <E>  and  vd  using  a  non-uniform  grid  were  compared  to  calculations  using  a  uni¬ 
form  grid  for  the  case  of  no  inelastic  processes  and  vM  T  =constant.  Since  analytic  solutions  to 
the  Boltzmann  equation  exist  for  the  case  of  vM  T  =  const,  these  parameters  could  then  be 
compared  to  the  expected  analytic  values.  Secondly,  vd,  and  a/N  were  calculated  on  both 
the  uniform  and  non-uniform  energy  grids  and  compared  to  experimental  data  for  several  gases, 
both  molecular  and  rare.  Finally,  with  the  electric  field  turned  off,  the  new  finite  differenced 
representation  for  the  momentum  transfer  recoil  term  (Eq  (A-17a)  and  Eq  (A- 17b))  was  tested 
for  convergence  to  the  proper  value  of  <E>  for  several  gas  temperatures. 

In  the  validation  of  the  non-uniform  energy  grid,  some  method  must  be  used  whereby  the  grid 
is  allocated.  For  this  analysis,  the  non-uniform  grid  was  established  by  using  a  pseudo- 
logaritinitic  method.  This  method  defines  the  width  of  each  energy  bin  by  the  relation 

Ae*  =  LnOfcr+E*.  (63) 

where  the  variables  E^,,  and  a  were  defined  so  that  the  maximum  energy  of  the  non-uniform 
energy  grid  was  approximately  the  same  as  the  maximum  energy  of  the  uniform  grid.  Thus,  the 
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energy  at  each  bin  of  the  non-uniform  grid  was  established  by  summing  the  width  of  each  bin 
below  it  When  any  comparisons  were  made  between  transport  parameters  calculated  using  the 
uniform  and  variable  energy  grids,  the  number  of  energy  bins  used  in  each  method  was  constant. 

Transport  Parameters  for  Constant  Collision  Frequency 

Since  an  analytic  solution  to  the  Boltzmann  equation  exists  for  the  special  case  of  a  constant 
momentum  transfer  collision  frequency,  this  condition  was  used  in  order  to  compare  the  calcu¬ 
lated  average  energy  and  drift  velocity  to  the  analytically  determined  ones.  The  analytic  values 
of  <E>  and  vd  were  determined  by  using  Eq  (44)  with  \MT  =  3.05  x  10'V  and  no  inelastic 
cross  sections  in  order  to  derive  the  reduced  electron  distribution  function,  /0(£).  From  this  func¬ 
tion,  the  average  energy  can  be  determined  by  Eq  (15).  The  drift  velocity  is  determined  by  using 
/0(e)  in  Eq  (44b)  to  get  /i(£),  from  which  vd  is  calculated  by  Eq  (13).  The  analytic  expressions 
derived  as  a  result  are 


<E>  =  4.1318  x10"2(E/N)2  eV  (64a) 

vd  =  1.412  x  105(E/N)  cm/s  (64'b) 

Table  1  shows  the  mparisons  between  <E>  calculated  with  a  uniform  and  variable  mesh  to 
the  analytic  values.  The  energy  axis  was  divided  into  150  bins,  with  the  maximum  energy  vary¬ 
ing  from  2  eV  at  2  Td,  to  200  eV  at  20  Td.  For  the  variable  energy  mesh,  Emin  and  a  were 
varied  so  that  the  maximum  energy  was  approximately  the  same  as  in  the  uniform  case.  The 
maximum  energy  for  the  calculation  was  determined  by  requiring  at  least  a  six  order  of  magni¬ 
tude  drop  in  the  distribution  at  the  maximum  energy,  but  not  more  than  a  ten  order  of  magnitude 
decrease.  Limiting  the  drop  to  less  than  ten  orders  of  magnitude  will  prevent  any  numerical 
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errors  arising  from  computer  roundoff.  Requiring  at  least  a  six  order  of  magnitude  drop  will 
ensure  that  nearly  all  of  the  electrons  are  accounted  for  in  the  calculation  of  the  transport  parame¬ 
ters. 


Table  1 

<E>  (eV)  for  vM  T  =  Constant,  Calculated  for  Uniform  vs  Non-Uniform  Energy  Grid 


E/N  (Td) 

Analytic 

Uniform 

%  Error 

Non-Uniform 

%  Error 

2 

0.165 

0.166 

0.606 

0.167 

1.212 

4 

0.662 

0.665 

0.453 

0.667 

0.755 

6 

1.487 

1.489 

0.135 

1.494 

0.471 

8 

2.644 

2.651 

0.265 

2.658 

0.530 

10 

4.132 

4.155 

0.557 

4.168 

0.871 

12 

5.950 

5.983 

0.555 

6.003 

0.891 

14 

8.098 

8.139 

0.506 

8.174 

0.939 

16 

10.577 

10.630 

0.501 

10.680 

0.974 

18 

13.387 

13.460 

0.545 

13.510 

0.919 

20 

16.527 

16.610 

0.502 

16.680 

0.926 

In  every  case,  the  code  has  over-estimated  the  average  energy.  The  values  calculated  with  the 
non-uniform  axis  are  consistently  higher  than  those  calculated  with  a  uniform  grid.  However,  it 
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must  be  remembered  that  the  variable  energy  grid  used  in  these  calculations  has  been  somewhat 
arbitrarily  constructed;  that  is,  each  energy  bin  width  was  assigned  according  to  Eq  (63).  No 
attempt  was  made  to  optimize  the  grid  for  the  calculations  of  drift  velcity  or  average  energy. 

Table  2  shows  a  similar  comparison  between  the  drift  velocity  calculations  with  the  different 
energy  grids  compared  to  the  analytic  values  for  different  E/N  values.  The  calculations  at  each 
E/N  were  accomplished  using  the  same  grid  mesh  that  was  used  previously  in  the  <E>  calcula¬ 
tions. 


Table  2 

Vd  (xlO6  cm/s)  for  vM  T  =  Constant,  Calculated  for  Uniform  vs  Non-Uniform  Energy  Grid 


E/N  (Td) 

Analytic 

Uniform 

%  Error 

Non-Uniform 

%  Error 

2 

0.2824 

0.2817 

0.248 

0.2813 

0.390 

4 

0.5648 

0.5633 

0.266 

0.5617 

0.549 

6 

0.8472 

0.8413 

0.696 

0.8413 

0.696 

8 

1.1296 

1.124 

0.496 

1.123 

0.584 

10 

1.4120 

1.409 

0.212 

1.406 

0.425 

12 

1.6944 

1.691 

0.200 

1.687 

0.437 

14 

1.9768 

1.971 

0.293 

1.966 

0.546 

16 

2.2592 

2.253 

0.274 

2.248 

0.496 

18 

2.5416 

2.535 

0.297 

2.529 

0.496 
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The  calculated  results  were  consistently  under-estimating  the  analytic  values,  with  the  uni¬ 
form  grid  generally  giving  better  estimations  of  the  velocities  than  the  non-uniform  grid.  Again, 
no  effort  was  made  to  construct  an  optimized  energy  grid  is  this  case. 

Figure  3  shows  a  comparison  between  the  energy  distributions  calculated  for  the  uniform  and 
non-uniform  grids  and  the  analytic  distribution  for  an  E/N  of  4  Td.  The  electron  distributions 
have  been  plotted  as  a  reduced  distribution  (cm“3eV~3/2)  such  that  a  Maxwellian  will  appear  as  a 
straight  line  when  plotted  on  a  semi-logarithmic  scale.  Plotted  on  this  scale,  there  is  very  little 
difference  between  each  of  the  distributions.  The  dotted  line  represents  the  distribution  calcu¬ 
lated  on  the  non-uniform  grid,  while  the  solid  line  represents  the  uniform  grid  calculation.  The 
dashed  line  represents  the  analytic  distribution  calculated  under  the  identical  conditions.  Both  of 
the  calculated  distributions  over-estimate  the  analytic  distribution  slightly,  illustrating  the  over¬ 
estimation  of  the  average  energy. 


Figure  3.  He  EDF  at  4  Td:  Uniform  vs  Non-Uniform 
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Figure  4  shows  the  percent  error  in  the  two  calculated  distributions  as  a  function  of  energy. 
This  error  is  defined  as  (/„, Jytic  VnumencX/uuiytk x  100. 

An  over-estimation  of  the  correct  number  density  will  result  in  a  negative  error,  while  an  under¬ 
estimation  will  result  in  a  positive  error.  At  low  energies,  the  non-uniform  grid  results  in  over¬ 
estimating  the  correct  values,  while  the  uniform  grid  under-estimates  the  correct  values.  At 
larger  energies,  both  of  the  calculated  distributions  over-estimate  the  analytic  result.  In  this 
region,  the  non-uniform  grid  over-estimates  the  analytic  distribution  more  than  the  uniform  grid, 
leading  to  a  greater  calculated  average  energy  for  the  former  compared  to  the  latter. 


Figure  4.  Error  in  Numerical  Calculation:  Uniform  vs  Non-Uniform 

The  finite  differenced  term  for  the  momentum  transfer  recoil  was  tested  for  convergence  to 
the  proper  functional  form.  In  order  to  do  this,  the  electric  field  was  zeroed,  and  the  momentum 
transfer  collision  frequency  was  held  constant  at  the  same  value  as  »hat  which  was  used  pre¬ 
viously.  The  maximum  energy  of  the  calculation  was  varied  for  each  gas  temperature  such  that 
the  distribution  underwent  at  least  a  six,  and  not  more  that  a  ten  order  of  magnitude  drop  between 
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the  lowest  and  highest  energies.  The  number  of  energy  bins  along  the  energy  axis  was  held  con¬ 
stant  at  200.  With  no  external  forces  applied  to  the  system  of  electrons,  the  distribution  should 
relax  to  a  Maxwellian  whose  average  energy  is  3/2  the  gas  temperature.  Table  3  illustrates  the 
results  for  a  variable  mesh  grid. 


Table  3 

<E>  (eV):  0  Electric  Field,  vM  r  =  Const,  Calculated  for  Non-Uniform  Energy  Grid 


Gas  Temp  (K) 

Analytic 

Calculated 

%  Error 

300 

3.876E-2 

3.819E-2 

1.471 

600 

7.752E-2 

7.780E-2 

0.361 

900 

1.163E-1 

1.168E-1 

0.430 

1200 

1.550E-1 

1.560E-1 

0.645 

1500 

1.938E-1 

1.948E-1 

0.516 

Figure  5  shows  the  converged  solution  plotted  as  a  reduced  distribution  for  a  Tg«  of  300  K.  In 
this  case,  a  was  0.0005  and  Emin  was  0.01  (Eq  (63)). 
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Energy  (eV) 

Figure  5.  EDF  Calculated  on  Non-Uniform  Energy  Grid  for  Tg„  =  300  K,  E/N  =  0 

Transport  Parameters  for  Gases 

Electron  EDF’s  were  calculated  for  Ar  and  H2  using  both  a  variable  and  uniform  energy  grid. 
The  resulting  transport  parameters  were  calculated  in  each  case  and  compared  to  experimental 
data  as  reported  by  Dutton  (5).  In  the  case  of  Ar,  the  parameters  were  drift  velocity,  average 
energy,  and  ct/N.  The  processes  considered  for  the  Ar  calculations  were  momentum  transfer, 
ground  state  ionization,  and  three  excited  states,  each  from  the  ground  state.  The  cross  sections 
for  these  processes  are  shown  in  Figure  6. 
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Figure  6.  Argon  Cross  Sections 

Figure  7  compares  the  drift  velocities  to  data  gathered  by  Dutton  (5).  Although  there  is  signifi¬ 
cant  scatter  in  the  data,  the  calculated  velocities  follow  the  general  trends  in  the  data.  The  data 
against  which  the  calculated  velocities  are  compared  was  measured  by  five  different  researchers, 
with  two  of  them  causing  the  largest  scatter  here. 
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4.50E+06' 


Figure  7.  Ar  Drift  Velocity  vs  E/N  (Data  from  Ref  5:612) 

Figure  8  shows  a  comparison  between  the  calculated  characteristic  energy  and  the  data,  which 
was  also  given  as  characteristic  energy.  Again,  the  general  agreement  is  fairly  good.  A  good 
agreement  between  the  calculated  and  experimental  values  of  energy  and  \d  indicate  that  the 
bulk  of  the  electron  distribution  has  been  accurately  computed.  This  is  due  to  the  characteristic 
dependence  of  both  these  parameters  on  the  bulk  of  the  electrons,  or  those  electrons  in  the  distri¬ 
bution  that  are  below  the  excitation  and  ionization  threshold  energies  of  the  gas. 
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Figure  8.  Ar  e k  vs  E/N  (Data  from  5:651) 

Cx/N  is  a  parameter  which  is  characteristic  of  the  tail  of  the  electron  distribution.  Thus,  if  the 

calculated  a/N  follows  the  experimentally  determined  values,  then  confidence  in  that  part  of  the 
distribution  governed  by  the  ionization  process  is  established.  Figure  9  shows  the  calculated  val¬ 
ues  of  a/N  and  the  experimental  data.  The  calculated  values  were  computed  without  any  super¬ 
elastic  process,  thus  the  slight  under-estimation  could  be  corrected  by  assuming  some  excited 
state  population  ratio  in  the  computation.  This  would  allow  the  superelastic  processes  to  drive 
electrons  out  to  the  higher  energies,  thus  increasing  the  ionization  rate  and  thereby  the  a/N.  This 
under-estimation  is  also  compounded  by  the  numeric  ionization  rate  calculation  method,  which 
in  the  code  is  a  simple  box  integrator  that  steps  backward.  This  backward  stepping  would  under¬ 
estimate  the  integration  of  the  distribution  with  the  cross  section,  thus  under-estimating  the  ion¬ 
ization  iate,  and  therefore  the  a/N.  This  effect  is  especially  prominent  if  the  distribution  falls  off 
rapidly  with  energy. 
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Figure  9.  Ar  aJN  vs  E/N  (Data  from  5:724) 

The  transport  parameters  calculated  in  H2  were  drift  velocity  and  Figures  10  and  1 1  show 

the  cross  sections  used  in  the  calculation  of  the  transport  parameters  for  molecular  Hydrogen. 
Figure  10  illustrates  the  momentum  transfer,  and  ground  state  electronic  and  ionization  cross  sec¬ 
tions.  Figure  1 1  shows  the  five  ground  state  vibrational  processes.  No  superelastic  processes 
were  used.  The  thrc  hold  energies  for  the  inelastic  processes  in  this  molecular  gas  are  in  sharp 
contrast  to  those  in  the  rare  gases,  such  as  Ar  or  He. 
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Figure  10.  H2  Cross  Sections  for  M.T.,  Excitation,  Ionization 


Figure  11.  H2  Cross  Sections  for  Ground  State  Vibrational  Excitation 
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The  drift  velocities  shown  in  Figure  12  follow  the  experimental  data  fairly  well  until  approxi¬ 
mately  30  Td,  beyond  which  the  calculated  values  are  greater  than  the  experimental.  The  veloci¬ 
ties  computed  with  the  uniform  and  non-uniform  energy  grids  are  in  close  agreement  throughout 
the  E/N  range  investigated,  indicating  that  this  error  does  not  lie  in  the  choice  of  energy  grids 
used. 


Figure  12.  H2  Drift  Velocity  vs  E/N  (Data  from  5:618) 


The  characteristic  energy  is  plotted  in  Figure  13.  In  this  plot,  the  calculated  average  energy 
has  been  multiplied  by  2/3  in  order  to  obtain  an  equal  weighting  with  the  data,  which  is  given  as 
characteristic  energy. 
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Figure  13.  H2  E*.  vs  E/N  (Data  from  5:655) 

In  comparing  the  transport  parameters  computed  with  the  uniform  and  variable  energy  grids 
for  H2  and  Argon,  it  appears  that  the  computed  data  for  the  molecular  gas  is  in  better  agreement 
with  the  data  than  the  rare  gas.  The  calculated  parameters  for  the  rare  gas  generally  diverges 
much  sooner.  This  result  is  probably  due  to  the  higher  inelastic  thresholds  in  Argon,  at  which 
point  the  logarithmic  grid  is  becoming  excessively  large  in  accurately  calculating  the  distribu¬ 
tion. 


Validation  of  Attachment  Process 


The  process  of  electron  attachment  was  validated  by  considering  the  effect  that  the  process 
would  have  on  the  electron  distribution,  as  well  as  verifying  the  predictions  made  by  Eq  (27)  for 
several  special  cases.  In  every  calculation,  energy  balance  was  considered  to  be  the  prime  com- 
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putational  diagnostic,  with  the  energy  always  balancing  to  within  0.1%.  In  this  validation,  the 
calculations  were  made  on  a  uniform  energy  grid  of  150  bins,  with  a  maximum  energy  of  20  eV. 
Argon  was  used  as  the  host  gas,  wherein  an  attaching  cross  section  was  added  of  varying  magni¬ 
tude  and  energy  dependence,  depending  on  the  case.  The  same  E/N  value  was  used  for  each  cal¬ 
culation.  The  other  processes  considered  were  momentum  transfer,  excitation,  and  ionization. 
Figure  14  shows  the  reduced  distriburion  resulting  from  the  inclusion  of  an  attachment  cross 
section  compared  to  the  no-attachment  case.  For  this  case,  the  attachment  cross  section  was  con¬ 
structed  such  that  the  attachment  rate  was  constant,  given  as  =  4.0  x  10"ncmV1.  The 

rate  was  defined  as  =  o(e)v(e)  fcmV]  at  each  energy.  With  a  constant  attachment 
frequency,  the  rate  at  which  electrons  will  be  removed  from  each  energy  bin  will  be  constant, 
thus  the  distribution  will  assume  the  same  form  as  if  attachment  were  not  considered.  In  as  much 
as  the  total  electron  number  density  is  smaller  in  the  case  in  which  attachment  is  considered,  a 
similar  form  translates  into  nxln2  =  constant.  As  expected,  the  form  of  the  distribution  remains 
similar,  although  the  distribution  in  which  attachment  was  considered  is  several  orders  of  magni¬ 
tude  lower  than  for  the  case  of  no  attachment.  The  calculated  average  energy  was  the  same  in 
each  case  (3.285  eV),  as  Eq  (27)  suggests. 
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Figure  14.  Calculated  EDF  with  a  Constant  Attachment  Rate 


Next,  an  attaching  cross  section  was  included  such  that  the  rate  of  attachment  increased  with 
energy.  The  actual  rate  used  rose  linearly  from  0.0  at  5  eV  to  2.0  x  10_®cm3s_1  at  20  eV.  This 
rate  was  given  as  =  l.OxlO^e  fore>5  eV  [cmV]. 

This  cross  section  results  in  electrons  at  the  higher  energies  having  a  greater  probability  of  being 
attached,  and  as  Figure  15  shows,  these  are  the  electrons  which  have  been  preferentially  removed 
from  the  distribution.  The  average  energy  for  this  distribution  was  calculated  to  be  1.792  eV, 
clearly  indicating  the  loss  of  the  higher  energy  electrons. 


61 


Figure  15.  Calculated  EDF  with  an  Attachment  Rate  Increasing  with  Energy 


Finally,  an  attaching  cross  section  that  decreased  rapidly  with  energy  was  used.  This  cross 
sec. .on  resulted  in  an  attaching  rate  that  also  decreased  with  energy,  from  a  maximum  of 
4.0  <  10~"  cmV  at  0.01  eV  to  0.0  at  2  eV.  The  EDF  with  and  without  attachment  are  shown  in 
Fig.  re  16.  The  loss  of  electrons  at  the  low  energies  is  readily  apparent,  and  in  fact,  this  portion 
of  le  curve  would  have  a  negative  temperature  associated  with  it.  This  form  of  the  attaching 
rate  will  shift  the  bulk  of  the  electrons  to  higher  energies,  resulting  in  an  average  energy  (5.134 
eV  that  is  greater  than  that  calculated  without  attachment  (3.285  eV). 
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Figure  16.  Calculated  EDF  with  an  Attachment  Rate  Decreasing  with  Energy 


Methods  of  Allocating  Variable  Energy  Grids 


Once  the  code  has  been  validated  with  the  use  of  a  non-uniform  energy  grid,  the  natural  q"es- 
tion  that  arises  is  "what  is  the  best  grid  to  use?".  A  logarithmically  allocated  bin  size  is  perhaps 
not  the  best  choice  for  Argon,  a  rare  gas  with  little  energy  separation  between  the  inelastic 
processes.  This  is  evidenced  by  the  disparity  between  the  transport  parameters  calculated  with 
uniform  and  non-uniform  energy  axes.  However,  the  use  of  the  log  grid  with  H2  appears  to  be 
practical  and  results  in  transport  calculations  that  are  similar  to  those  calculated  with  the  uniform 
grid.  It  would  be  advantageous  if  some  method  could  be  devised  in  which  the  non-uniform 
energy  axis  is  defined  so  it  would  increase  the  accuracy  of  the  calculation  of  certain  parameters 
over  those  calculated  with  a  uniform  grid,  for  the  same  number  of  grid  points.  These  parameters 
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could  include  normalization  of  the  distribution,  drift  velocity,  average  energy,  and  kinetic  rate 
calculation.  The  first  three  of  these  parameters  depends  predominately  on  the  bulk  of  the  distri¬ 
bution,  while  the  last  depends  mostly  on  the  tail  of  the  distribution.  Thus,  while  a  non-uniform 
grid  that  concentrates  bins  in  the  bulk  may  increase  the  accuracy  of  the  drift  velocity,  average 
energy,  and  normalization  integrals  over  those  calculated  with  a  uniform  grid,  the  resulting 
kinetic  rate  information  may  be  much  worse  than  the  uniform  calculation.  Likewise,  a  grid 
which  places  many  bins  in  the  tail  of  the  distribution  may  result  in  a  rate  calculation  that  is  an 
improvement  over  the  uniform  grid,  but  now  the  other  three  parameters  are  far  off  the  mark.  The 
ideal  case  would  be  a  grid  allocated  .uch  that  the  accuracy  of  all  parameters  is  improved.  This  is 
the  motivation  for  Eq  (42),  with  an  intelligent  choice  for  the  function  G. 

In  order  to  help  explore  this  possibility,  Eq  *44)  was  used  as  a  distribution  function.  The 
simplifying  assumption  that  f0  doesn’t  change  much  over  the  interval  e  to  e  +  5e*  allows  the 
imbedded  fa  terms  on  the  right  hand  side  of  Eq  (44)  to  cancel  one  another.  This  assumption, 
while  not  strictly  valid  for  regions  around  the  threshold  energy,  will  nevertheless  allow  Eq  (44) 
to  serve  as  a  starting  point  in  the  dynamic  mesh  allocation  method.  With  the  previous  assump¬ 
tion,  Eq  (44)  reduces  to  a  more  tractable  Eq  (65). 

e 

/»<£)  =  eXP['(E^)j(^/XQ"(X)Q'(X)dX 

v  '  o 

c  ^  ,  e  +  &a 

rQ/(x)  c 

+  J—  J  Q>,(y)ydydx)]  (65) 

0  E 

where  f„  is  the  reduced  distribution.  Q*  will  be  assumed  to  be  a  step  function,  with  some  thresh¬ 
old  energy  below  which  the  cross  section  is  zero,  and  above  which  the  cross  section  is  con- 
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stant  The  right  hand  side  of  Eq  (65)  is  then  analytic,  and  can  be  used  to  represent  the 
approximate  solution  to  the  Boltzmann  equation,  with  which  the  validating  parameters  are 
calculated. 


e 

/”<e)  =  exp["(i^j5(^  JxQ"(x)a(x)dx 

+  J-^— J  Q*(y)ydydx)]  (66) 

0  E 

The  parameters  used  were  normalization  of  the  distribution,  average  energy,  drift  velocity  and 
excitation  rate.  Although  the  distribution  function  itself  was  analytic,  the  integrations  required 
for  these  comparison  parameters  could  not  be  carried  out  analytically  when  a  non-zero  inelastic 
cross  section  was  included  in  Eq  (65).  Therefore,  these  integrations  were  accomplished  numer¬ 
ically,  using  a  uniform  energy  grid  of  400  points,  and  used  as  the  "standard".  The  other  uniform 
and  non-uniform  integrations  were  made  on  a  50  point  grid,  and  compared  to  this  standard. 

Several  choices  for  the  function  G  in  Eq  (42)  were  tried,  with  varying  degrees  of  success  in 
improving  the  calculation  of  the  normalization,  drift  velocity,  and  kinetic  rate  over  the  uniform 
integration.  The  choices  for  these  merit  functions  were  driven  by  their  respective  influence  on 
the  placement  of  energy  bins  on  the  dynamic  energy  axis.  Since  an  allocation  function  which 
would  place  more  bins  in  the  body  of  the  distribution  wouldn’t  aid  in  the  rate  calculation,  and 
vice-versa,  two  merit  functions  were  used  to  work  in  conjunction  with  each  other.  In  this 
method,  the  final  merit  function  was  the  average  of  two  others,  one  designed  to  aid  the  calcula¬ 
tion  of  bulk  dependent  parameters  (average  energy,  drift  velocity)  and  the  other  designed  to  aid 
in  the  tail  dependent  ones  (kinetic  rates).  The  function  chosen  for  the  bulk  optimization  was  an 
energy  moment  of  the  distribution,  where  G  was  represented  by 
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G,(e)  =  A“/.(e)  (67) 

where  n  is  some  power.  If  n  =  0,  then  G,  would  be  a  normalization  integrand.  If  n  =  1,  G,  would 
be  an  average  energy  integrand.  Increasing  the  power  of  n  will  shift  the  peak  of  G,  out  to  higher 
energies,  which  will  delay  the  peak  of  the  cumulative  integral  function  F.  In  testing  this  algo¬ 
rithm,  it  was  found  that  the  value  of  n  that  would  result  in  the  best  integration  of  the  average 
energy  and  drift  velocity  for  the  non-uniform  grid  was  dependent  on  the  E/N  ratio. 

The  integration  of  the  tail  was  aided  by  constructing  a  function  G  given  by  the  rate  integrand: 

G2(e)  =  Q(e)e/0(e)  (68) 

Thus  G  =  Gj  +  G2. 

In  integrating  the  rate  parameter,  it  became  apparent  that  the  exact  placement  of  the  threshold 
energy  was  a  very  critical  variable.  Consider  a  uniform  energy  grid  where  the  threshold  energy 
is  just  slightly  less  than  the  bin  energy,  as  in  Figure  17.  In  this  instance,  the  bin  energy  is  defined 
as  that  energy  corresponding  to  the  right  hand  side  of  a  particular  finite  differenced  bin.  If  the 
inelastic  cross  section  rises  steeply,  and  the  distribution  is  falling  off  rapidly,  then  the  resulting 
rate  integrand  will  appear  as  that  shown.  The  simple  Simpson  backward  integrator  would  result 
in  a  large  over-estimation  of  the  integral  from  the  first  bin  containing  a  non-zero  value  for  the 
rate.  On  the  other  hand,  the  backward  integration  of  the  other  bins  would  always  result  in  the 
under-estimation  of  the  actual  value  of  the  integral.  The  result  of  these  two  effects  is  that  the 
calculated  integral  may  be  high,  low,  or  very  close  to  the  actual  value  of  the  integral.  If  the  over¬ 
estimation  is  approximately  the  same  as  the  under-estimation,  then  the  calculated  integration  will 
be  close  to  the  correct  value.  This  result  will  occur  in  spite  of  the  errors  in  the  integration 
method. 
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Figure  17.  Rate  integrand  with  slightly  less  than  a  bin  energy 

Consider  now  a  threshold  energy  slightly  greater  than  a  bin  energy,  as  in  Figure  18.  In  this 
case,  the  rate  integral  is  only  under-estimated,  and  could  be  worse  than  the  value  calculated  by 
the  previous  method. 

With  the  non-uniform  energy  grid  construction  technique  used  here,  the  bin  before  the  thresh¬ 
old  energy  is  typically  fairly  wide,  resulting  in  a  significant  over-estimation  of  the  rate  integral. 
This  is  especially  so  with  the  addition  of  many  more  bins  into  the  region  due  to  ♦h?  more 
accurate  integration  in  the  following  bins.  Therefore,  the  placement  of  a  bin  at  the  threshold 
energy  becomes  an  important  step  in  accurate  rate  calculation. 
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Figure  18.  Rate  integrand  with  slightly  greater  than  a  bin  energy  on  uniform  energy  grid 

Figure  19  shows  the  result  of  using  this  choice  for  an  energy  bin  allocation  function  in  defin¬ 
ing  a  50  point  non-uniform  grid.  In  this  plot,  the  normalized  distribution  function  and  normal¬ 
ized  rate  integrand  are  shown  in  order  to  visualize  not  only  the  construction  of  the  allocation 
function  itself,  but  also  in  order  to  understand  the  resulting  integrations  of  each  over  the  variable 
energy  axis,  which  is  shown  as  the  variable  size  grid.  The  inelastic  threshold  is  at  an  energy  of 
5.06. 
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Figure  19.  EDF  and  Rate  Integrand  and  resulting  Allocation  Function 

It  is  instructive  to  see  how  the  variable  energy  grid  will  integrate  each  of  the  parameters  of 
interest.  Figure  20  shows  each  of  these  parameters  in  a  box  layout.  At  the  top  left  is  shown  the 
energy  moment  ot  the  distribution  along  with  the  energy  grid  resulting  from  the  merit  function  of 
G.  The  integral  of  this  function  yields  the  average  energy.  Below  the  energy  moment  is  the  rate 
integrand.  The  increase  of  the  bin  density  in  this  region  always  resulted  in  a  more  accurate  rate 
calculation  compared  to  the  uniform  method.  The  top  right  panel  illustrates  how  the  variable 
axis  will  integrate  the  anisotropic  part  of  the  distribution,  which  yields  the  drift  velocity.  The 
bottom  right  plot  shows  the  symmetric  part  of  the  distribution  function,  whose  integration  will 
result  in  the  normalization  that  is  used  in  both  the  rate  and  drift  velocity  calculations. 
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Figure  20.  Non-Uniform  Energy  Grid  with  Validation  Parameters 


Figure  21  shows  the  allocation  function  G,  where  n  =  1.0  in  Eq  (66),  at  5  Td  and  the  distribu¬ 
tion  and  rate  functions  from  which  it  is  constructed.  Comparing  this  plot  to  Figure  22,  which  is 
for  the  same  conditions  except  n  =  0.2,  leads  to  some  interesting  points.  As  the  value  of  n 
increases,  its  effect  on  the  allocation  function  is  evident  where  the  slope  is  reduced  at  the  lower 
energies,  resulting  in  a  less  dense  allocation  of  energy  bins  in  this  region.  Varying  the  value  of  n 
has  little  effect  on  the  allocation  of  energy  bins  at  the  higher  energies  where  the  inelastic  cross 
section  lies,  thus  the  calculation  of  the  rate  varies  little. 
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Figure  21.  Allocation  Function  for  5  Td,  n  =  1.0 
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Figure  23  sho-’  s  a  comparison  of  the  validation  parameters  for  different  values  of  n,  with  the 
uniform  calculations  for  the  same  number  of  bins.  It  was  found  that  the  optimum  value  of  n  in 
Eq  (66)  for  5  Td  was  0.6.  At  n  values  lower  than  0.6,  the  calculated  normalization,  average 
energy  and  drift  velocity  parameters  were  all  under-estimating  the  standard,  while  higher  values 
led  to  an  over-estimation.  Varying  the  n  value  had  a  very  small  effect  on  the  non-uniform  rate 
calculation,  which  was  always  better  than  the  uniform  calculation,  and  is  not  included  in  this 
plot. 


72 


As  the  E/N  ratio  is  reduced,  the  distribution  becomes  more  localized  toward  the  lower  ener¬ 
gies,  which  results  in  the  optimum  n  value  increasing,  as  seen  in  Figure  24.  Figures  25  and  26 
illustrate  how  the  different  n  values  lead  to  their  respective  integrations  of  the  validation 
parameters.  At  the  higher  n  value  (Figure  25)  the  smaller  slope  of  the  allocation  function  at  low 
energies  allows  more  bins  to  placed  in  the  region  where  the  distribution  is  falling  off  rapidly. 
Although  this  movement  of  energy  bins  to  the  higher  energies  results  in  a  poorer  integration  of 
the  low  energy  portion  of  the  distribution,  it  is  offset  by  the  increased  accuracy  at  the  higher 
energies.  The  lower  n  value  (Figure  26)  leads  to  too  many  points  allocated  near  the  origin,  thus 
under-estimating  the  integration  of  the  tail. 
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Figure  24.  n  vs  Percent  Error  for  4  Td 
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Figure  25.  Allocation  Function  for  4  Td,  n  =  1.0 
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Figure  26.  Allocation  Function  for  4  Td,  n  =  0.2 

It  is  interesting  to  note  that  even  with  the  optimum  value  of  n  used  for  the  non-uniform  calcu¬ 
lations,  a  3.8%  increase  in  the  accuracy  of  the  rate  comes  at  the  price  of  a  38.8%  increase  in  the 
number  of  bins  in  the  cross  section  region. 

Electron-Electron  Collisions  with  Non-Uniform  Energy  Grid 

It  is  well  known  that  electron-electron  collisions  tend  to  Maxwellianize  the  distribution,  thus 
smoothing  out  any  sharp  features  that  tend  to  result  from  inelastic  collisions.  In  the  absence  of 
any  other  effects,  electrons  colliding  with  themselves  will  relax  toward  a  Maxwellian  with  an 
electron  temperature  equal  to  2/3  of  the  average  energy.  The  electrons  must  obey  this  property 
regardless  of  the  initial  distribution.  Conservation  of  energy  and  of  particles  must  also  be  obeyed 
with  each  iteration. 

The  electron-electron  interactions  on  the  non-uniform  energy  grid  were  tested  in  several 
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ways.  With  an  initial  distribution  of  a  Maxwellian  with  a  temperature  of  2/3  the  average  energy, 
the  iterative  solution  was  allowed  to  propagate  forward  in  time.  The  test  for  the  distribution  con¬ 
vergence  was  disabled,  allowing  each  time  iteration  to  proceed  until  a  fixed  number  of  loops 
through  the  electron  collision  subroutine  were  accomplished,  in  this  case  2000.  Figure  27  shows 
the  resulting  distribution,  along  with  the  initial  Maxwellian  distribution. 


Figure  27.  EDF  with  electron-electron  collisions  only 

The  initial  distribution  was  observed  to  depart  from  the  Maxwellian  as  the  program  iterated 
through  the  solutions.  This  departure  is  more  evident  at  the  lower  energies.  Conservation  of 
energy  and  of  particles  were  obeyed  at  each  iteration,  within  the  limits  of  the  numerical  represen¬ 
tation  of  each.  When  this  test  case  was  performed  on  the  uniform  grid,  the  initial  and  final  distri¬ 
bution  were  exactly  the  same  out  to  at  least  six  decimal  places,  with  conservation  of  energy  and 
particles  observed.  The  deviation  from  the  Maxwellian  with  the  algorithm  using  a  non-uniform 
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grid  indicates  that  the  At ,  matrix  is  not  properly  constructed.  Since  the  particles  are  conserved, 
the  boundary  conditions  applied  to  this  matrix  seem  to  be  correct.  The  relationship  between  the 
AkJ  and  Bt ,  matrices  also  appears  to  be  correct  due  to  energy  conservation. 

In  order  to  test  the  algorithm's  ability  to  relax  to  the  Maxwellian  form,  an  initial  distribution 
was  constructed  with  3.6  x  1012  electrons  at  1.5  eV,  and  a  small  fraction  (1.0  x  105)  in  every  other 
bin.  The  temperature  of  the  resultant  Maxwellian  must  then  be  very  close  to  1  eV.  Figure  28 
shows  the  distribution  at  various  points  throughout  the  iteration  sequence.  Since  the  time  step 
used  was  0.2  ns,  the  distributions  at  500,  1000,  and  1500  iterations  correspond  to  times  of  100  ns, 
200  ns,  and  300  ns  respectively. 


Figure  28.  EDF  with  electron-electron  only;  spike  initial  distribution 

Although  it  appears  that  the  distribution  is  slowly  relaxing  to  the  Maxwellian  at  an  electron 
temperature  of  1  eV,  conservation  of  energy  and  of  particles  was  not  obeyed  to  all  orders.  The 
number  density  of  electrons  was  observed  to  increase  with  time,  causing  the  system  to  therefore 
gain  energy  with  time.  This  breakdown  indicates  a  probable  error  in  the  implementation  of  the 
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boundary  conditions,  even  though  conservation  of  particles  was  obeyed  in  the  previous  test  case. 
Even  though  the  relationship  of  the  AM’s  and  B4  ,’s  was  designed  such  that  conservation  of 
energy  is  always  obeyed,  this  cannot  hold  if  the  particles  are  not  conserved.  When  this  same  test 
case  was  performed  with  the  uniform  grid,  conservation  of  energy  and  particles  was  observed  at 
each  iteration,  with  the  distribution  approaching  the  1  eV  Maxwellian  in  approximately  the  same 
manner  as  Figure  28. 

Electron  Beams 


In  order  to  isolate  the  effect  of  the  electron  beam  alone,  the  electric  field  was  turned  off, 
thereby  eliminating  any  field  driven  flux.  A  logarithmic  energy  axis  was  used,  and  a  point 
source  of  electrons  was  placed  at  60.4  eV.  The  maximum  energy  of  the  grid  was  just  slightly 
greater  than  this  value.  A  beam  current  of  10  Amps  was  introduced  into  the  H2  gas,  which  was  at 
1  atmosphere  and  300  K.  The  processes  considered  were  vibrational  excitations  from  the  ground 
state  to  levels  1  through  5,  the  first  electronic  excitation  from  the  ground  state,  and  ground  state 
ionization.  Figure  29  shows  the  converged  solution  26  ns  after  the  beam  turned  on,  with  a  Max¬ 
wellian  being  the  initial  distribution.  For  the  electrons  being  injected  into  the  plasma  at  60  eV, 
there  are  only  two  places  to  go.  They  can  either  excite  the  ground  state  to  the  first  electronic 
state  (12.6  eV),  or  they  can  ionize  the  ground  state  H2  molecule  (15.427  eV).  The  large  increase 
in  population  at  47  eV  is  due  to  the  first  process,  while  the  peak  at  44.5  eV  is  due  to  the  second. 
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Figure  29.  H2  EDF  with  Electron  Beam  Source  at  60.4  eV 


In  the  calculation  of  this  distribution  function,  it  was  necessary  to  zero  all  momentum  transfer 
flux  out  of  the  energy  bin  containing  the  source  beam.  Without  doing  so  resulted  in  a  distribu¬ 
tion  centered  around  the  source  that  was  widely  spaced  in  energy,  covering  r.t»proximately  25  eV. 
This  broad  region  was  due  to  the  large  downflux  of  electrons  due  to  recoil  with  the  heavy 
particles.  Additionally,  in  the  region  between  the  secondary  peaks  and  the  source  it  was  neces¬ 
sary  to  set  a  lower  limit  on  the  number  density  of  electrons  in  the  energy  bins.  Due  to  the  con¬ 
stant  depletion  of  electrons  in  this  region  without  any  influx  of  electrons  from  either  below  or 
above,  the  population  will  continue  to  decrease  with  each  iteration  of  the  time -dependent 
solution.  As  a  result,  the  computation  went  unstable  (ie,  negative  number  densities)  when  the 
computer  reached  its  smallest  representable  number,  which  was  of  the  order  of  lO"300.  In  order  to 
avoid  this  from  happening,  an  articial  lower  limit  of  1  x  10'50  was  placed  on  the  electron  densities 
at  these  energies. 
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The  distribution  is  relatively  flat  across  the  intermediate  region  due  to  the  small  energy  cou¬ 
pling  between  the  H2  molecules  and  the  electrons,  which  is  2m/M.  In  this  region  momentum 
transfer  is  the  dominate  energy  transfer  mechanism.  The  electrons  will  continue  to  exchange 
energy  with  the  heavy  particles  in  this  manner  until  the  vibrational  cross  sections  become  effec¬ 
tive.  At  the  very  low  energies  (approximately  2  eV  and  below),  the  gas  molecules  at  300  K  are 
again  exchanging  energy  with  the  electrons,  with  the  distribution  becoming  a  Maxwellian  at 
appromimately  this  same  temperature. 

In  comparing  this  distribution  to  that  calculated  by  Bretagne  (3:821),  the  same  general  fea¬ 
tures  are  observed.  However,  there  is  much  more  structure  in  Bretagne’s  calculation  of  the  dis¬ 
tribution.  It  appears  that  the  first  spike  in  the  his  distribution  below  the  source  beam  (50  eV)  is 
due  to  electronic  excitation  in  H2,  either  Is  to  2p  or  Is  to  2s,  both  of  which  have  a  10.2  eV 
threshold.  The  large  population  just  below  this  (45  eV  or  so)  is  probably  due  to  ionization.  The 
gap  that  appears  between  these  two  peaks  is  likely  explained  by  Bretagne’s  consideration  of  elec¬ 
tron  loss  mechanisms,  specifically  recombination  and  diffusive  wall  losses.  These  effects  were 
not  considered  in  the  present  calculation,  and  the  resultant  lack  of  structure  between  the  two 
peaks  in  Figure  29  is  due  to  a  large  momentum  transfer  downflux  from  the  energy  bin  filled  by 
electronic  excitation.  Likewise,  any  structure  that  may  appear  in  the  intermediate  region  is 
"washed"  smoothly  away  due  to  this  same  effect.  By  introducing  loss  mechanisms,  not  all  of  the 
electrons  will  be  able  to  collide  elastically  with  there  collision  partners,  since  many  will  be  lost 
out  of  the  continuum  due  to  the  recombination  or  diffustion.  As  a  result,  the  downflux  of 
momentum  transfer  from  one  energy  bin  to  the  next  lower  bin  will  be  reduced,  and  the  structure 
that  Bretagne  reports  may  then  be  observable. 

The  calculated  solution  with  the  source  at  90.0  eV  is  shown  in  Figure  30.  The  basic  structure 
of  the  distribution  is  similar  to  that  with  the  60.4  eV  source.  Here  again,  the  differences  between 
Bretagne  and  the  calculated  distribution  lie  in  the  lack  of  detailed  structure. 
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Figure  30.  H2  EDF  with  Electron  Beam  Source  at  90.0  eV 


Negative  Total  Electron  Mobility 

In  order  to  explore  -  :  possibility  of  negative  mobility,  the  solution  method  utilizing  the  time- 
dependent  Boltzmann  equation  was  exercised  under  various  conditions  of  physical  interest. 

Since  the  principle  means  by  which  a  negative  mobility,  and  therefore  a  negative  drift  velocity, 
can  occur  is  a  low  energy  loss  of  electrons,  this  process  in  the  form  of  attachment  was  included. 
Fluorine,  which  has  a  large  attachment  cross  section  at  low  energies,  was  used  as  the  attaching 
gas,  with  Argon  used  as  the  primary  gas.  Several  gas  concentration  levels  were  explored,  result¬ 
ing  in  a  strong  influence  on  the  drift  velocity. 

The  momentum  transfer  and  attachment  cross  section  for  Fluorine  are  shown  in  Figure  31. 
The  large  cross  section  for  attachment  at  the  low  energies  results  in  a  high  probability  of  the 
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electrons  to  be  lost  at  these  same  energies.  Since  the  cross  section  falls  off  very  rapidly,  the  elec¬ 
trons  at  energies  greater  than  2  eV  or  so  will  be  largely  unaffected,  with  the  result  being  a 
reduced  distribution  with  a  positive  slope  at  low  energy,  and  a  more  traditional  negative  slope  at 
the  higher  energies. 


Figure  31.  Fluorine  Cross  sections 

The  energy  distribution  calculations  were  performed  on  a  150  point  uniform  energy  grid,  that 
was  typically  8  eV  at  the  maximum  energy.  Initially,  a  gas  mix  of  100%  Ar  and  0%  F2  was  used 
to  validate  the  drift  velocities  calculated  from  the  code.  The  calculated  drift  velocities  were  com¬ 
pared  with  data  (5:612)  resulting  in  fair  agreement  (Figure  32). 
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Figure  32.  Drift  Velocity  in  Argon  at  Low  E/N 

Two  mixes  of  Ar/F2  were  used,  the  first  being  99.9%  Ar  with  0  1%  F2,  and  the  second  °9.5% 
Ar  and  0.5%  F2.  At  low  E/N,  the  field  driven  flux  is  very  weak,  resulting  in  a  distribution  that 
falls  off  rapidly  at  higher  energies.  This  rapid  decrease  in  the  distribution  is  due  to  the  small 
field  driven  flux  combined  with  the  Argon  momentum  transfer  cross  section  that  increases  rap¬ 
idly  with  energy,  resulting  in  a  collision  frequency  that  increases  rapidly.  As  this  collision  fre¬ 
quency  increases,  the  elastic  energy  transfer  becomes  greater,  resulting  in  the  steep  negative 
slope.  At  the  low  energies,  the  attachment  process  will  continually  serve  as  an  electron  sink, 
resulting  in  a  distribution  with  low  populations  in  this  energy  region.  The  EDF  for  a  99.5/0.5 
mix  of  Ar/F2  at  0.5  Td  is  shown  in  Figure  33.  The  low  energy  region  with  the  positive  slope  is 
responsible  for  the  negative  drift  velocity.  This  region  also  has  a  characteristic  temperature  that 
is  negative. 
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The  resulting  drift  velocities  are  plotted  against  E/N  in  Figure  34.  The  distribution  with  the 
higher  concentration  of  Fluorine  requires  a  greater  upflux  of  electrons  in  order  to  recover  a  posi¬ 
tive  drift  velocity.  Since  this  upflux  is  caused  by  the  electric  field,  a  larger  E/N  ratio  is  required 
to  obtain  a  positive  drift  velocity  compared  to  the  lower  concentration  casey 
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Figure  34.  Drift  Velocity  in  Ar/F2  vs  E/N 

In  both  mixtures,  as  the  E/N  is  increased,  the  field  driven  flux  increases,  while  the  effect  of 
the  attachment  process  is  reduced.  This  results  in  the  calculated  drift  velocities  for  both  mixes 
approaching  that  of  pure  Argon,  with  the  lower  F  concentration  case  approaching  it  much  sooner. 
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VI.  Conclusion  and  Recommendations 


The  code  developed  by  Seger  (16)  has  been  improved  to  include  the  processes  of  superelas¬ 
tics,  attachment,  and  excitation  and  ionization  from  excited  states.  The  capability  to  solve  the 
Boltzmann  equation  on  a  non-uniform  energy  grid  has  also  been  added  to  the  code.  This  added 
capability  is  essential  in  order  to  extend  the  energy  range  of  the  solution  past  a  previous  upper 
limit  of  30  eV  or  so.  With  the  non-uniform  energy  grid,  the  upper  limit  has  been  increased  three¬ 
fold,  to  a  demonsuated  90  eV.  These  larger  energy  ranges  are  essential  in  order  to  explore  the 
effects  of  an  electron  beam  on  the  energy  distribution  functions,  which  have  been  calculated 
under  conditions  typical  for  a  hydrogen  magnetic  multicusp  discharge.  In  addition,  the  recoil 
term  for  momentum  transfer  has  been  finite  differenced  using  a  new  method  that  recovers  the 
expected  Maxwellian  form  in  the  absence  of  electric  fields,  as  well  as  maintaining  numerical  sta¬ 
bility  under  this  same  condition,  both  of  which  were  lacking  in  the  previous  version  of  the  code. 
Transport  parameters  computed  with  the  non-uniform  energy  grid  are  generally  in  agreement 
with  those  calculated  with  the  uniform  grid  as  well  as  experimental  data.  The  use  of  the  loga¬ 
rithmic  grid  yields  more  accurate  computations  of  the  transport  parameters  for  H2  than  for  Ar. 
This  is  likely  due  to  the  grouping  of  inelastic  processes  at  the  lower  energies  in  H2,  whereas  in  Ar 
the  inelastic  processes  occur  at  higher  energies.  It  is  expected  that  the  accuracy  of  the  calcula¬ 
tion  of  the  transport  parameters  with  a  log  grid  in  He  would  be  worse  than  in  Ar,  due  to  its 
greater  inelastic  thresholds.  The  optimization  of  the  construction  of  a  non-uniform  grid  has  been 
explored  with  limited  success.  A  more  accurate  kinetic  rate  calculation  is  quite  feasible, 
although  it  seems  that  the  price  for  such  comes  in  the  form  of  a  less  accurate  calculation  of  the 
drift  velocity,  average  energy,  and  normalization  -  parameters  which  depend  mainly  on  the  bulk 
of  the  distribution.  The  benefit  of  placing  an  energy  bin  at  the  threshold  value  of  each  inelastic 
process  is  established,  resulting  in  a  more  accurate  kinetic  rate  calculation.  The  utility  of  this 
"threshold  marking"  technique  is  enhanced  when  a  cross  section  with  a  large  slope  at  threshold  is 
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used.  Finally,  the  possibility  of  steady  state  negative  drift  velocities  has  been  explored  in  lean 
Ar/F2  mixtures.  The  numerical  calculations  indicate  that  at  low  E/N  values  and  with  the  proper 
gas  concentrations,  negative  drift  velocities  can  be  achieved. 

While  these  added  capabilities  increase  the  ability  to  analyze  plasmas,  more  can  and  must  be 
done  to  increase  the  accuracy  of  the  solution,  as  well  as  the  transport  parameters  computed  from 
it.  A  better  method  for  establishing  the  non-uniform  axis  is  essential  in  order  to  achieve  its 
potential  at  being  more  accurate  than  the  uniform  axis.  As  shown  previously,  the  logarithmic 
grid  is  better  in  molecular  gases  than  in  rare  gases,  with  the  uniform  axis  being  generally  better 
in  either.  However,  it  is  still  maintained  that  with  the  proper  grid  the  non-uniform  can  be  better 
than  uniform,  at  least  in  the  calculation  of  some  of  the  parameters  of  interest  The  optimization 
of  this  non-uniform  grid  needs  to  be  explored  more  fully  in  order  to  achieve  this  goal. 

In  electron  beam  generated  plasmas,  the  inclusion  of  recombination  and  diffusive  wall  losses 
appears  to  be  important  in  reproducing  the  structure  that  the  distributions  calculated  by  Bretagne 
(3)  exhibit.  Additionally,  electron  -  electron  interactions  within  the  electron  beam  generated 
plasma  will  tend  to  Maxwellianize  the  distribution,  and  thereby  serve  to  smooth  out  any  sharp 
structure  arising  in  the  distribution  from  inelastic  collisions.  Therefore,  its  inclusion  into  the 
solution  method  is  important,  especially  at  the  higher  fractional  ionization  levels. 

Another  problem  with  the  present  code  exists  in  the  finite  differenced  terms  for  the  field 
driven  flux  when  used  in  conjuction  with  an  e-  beam.  As  has  been  reported,  very  sharp  structure 
can  result  in  the  distribution  as  a  consequence  of  the  electrons  being  fed  at  high  energies,  and 
relaxing  to  lower  energies  through  inelastic  collisions.  The  field  driven  flux  term  includes  a  first 
and  second  derivative  of  the  electrons  number  density  with  respect  to  energy.  Since  any  sharp 
features  in  the  distribution  will  cause  both  of  these  derivatives  to  be  discontinuous,  numerical 
instabilities  around  these  regions  have  resulted  when  the  present  code  is  used  with  both  an  elec¬ 
tric  field  and  an  electron  beam.  It  may  be  possible  to  eliminate  part  or  all  of  this  instability  by 
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finite  differencing  the  field  driven  flux  term  to  higher  orders  and/or  using  a  central  difference 
scheme  in  place  of  a  forward  difference  one.  Doing  so,  however,  will  result  in  the  C  matrix  no 
longer  being  tridiagonal,  which  will  increase  the  time  required  to  decompose  it. 

In  order  to  isolate  the  electron  beam  generated  plasma  from  the  presence  of  an  electric  field, 
the  recoil  momentum  transfer  flux  term  had  to  be  finite  differenced  in  a  new  way,  based  on  that 
used  by  Greene  and  Elliott  (6).  With  no  electric  field,  and  in  the  absence  of  any  source,  this 
version  of  the  recoil  flux  term  has  maintained  numerical  stability  and  reproduced  a  Maxwellian 
at  a  temperature  of  2/3  of  the  average  gas  temperature  as  required  by  physical  laws.  However, 
when  this  representation  of  the  recoil  term  was  used  in  conjunction  with  an  electric  field  in 
Argon,  the  distribution  was  overpopulated  at  the  lower  energies,  and  underpopulated  at  the 
higher  ones.  The  resultant  transport  parameters  were  also  well  off  the  experimental  values. 
Returning  to  Rockwood’s  recoil  term  re-derived  for  the  non-uniform  energy  grid  (Appendix  A) 
resulted  in  the  proper  calculation  of  the  distribution  function  and  its  transport  parameters.  How¬ 
ever,  when  this  version  of  the  recoil  term  is  used  in  the  absence  of  an  electric  field,  numerical 
instability  is  the  result.  It  would  seem  likely  that  there  is  some  value  of  the  electric  field  at  which 
Elliott  and  Greene’s  version  will  produce  correct  results.  More  testing  and  investigation  is 
required  in  order  to  determine  where  the  cutoff  between  one  representation  and  the  other  will  be. 
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Appends.  A 


Derivation  of  field  and  recoil  flux  terms  for  a  non-uniform  energy  axis 

Although  the  terms  representing  the  flux  of  electrons  through  energy  space  due  to  the  applied 
field  and  recoil  have  been  derived  in  the  past,  the  inclusion  of  a  non-uniform  energy  axis  requires 
that  they  be  rederived  explicity.  In  this  section,  they  will  be  considered  separately,  and  then 
grouped  together  at  the  end.  First,  consideration  will  be  given  to  the  field  driven  flux  term,  with 
the  momentum  transfer  term  following.  The  energ'  axis  will  be  define  as 
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Figure  35.  Energy  Axis  for  Momentum  Transfer 

Field  Driven  Flux 

The  origin  of  the  field  driven  flux  term  as  presently  implemented  in  Megaboltz  stems  from  a 
paper  by  Zel’Dovich  and  Raizer  (19),  as  reported  by  Greene  (7).  In  this  paper,  they  expressed 
the  continuity  of  electrons  under  the  influence  of  an  intense  electromagnetic  pulse  as 
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where  J  is  given  by 


J  =  /l(£)u- 

where  n(e)  is  the  number  density  of  electrons  with  energy  between  e  and  e+  6e,  u  is  the  "veloc¬ 
ity"  of  the  election  along  the  1  dimensional  energy  axis,  and  D  is  the  diffusion  coefficient  along 
the  energy  axis.  The  velocity  of  the  electron  along  the  energy  axis  can  be  envisioned  as  being 
similar  to  the  velocity  of  the  electron  through  phase  space.  In  such  an  analogy,  the  geometric 
divergence  of  the  flux  («(v)v)  through  phase  space  is  equivalent  to  the  energy  divergence  of  the 
flux  (n(e)u)  through  energy  space.  Here,  the  velocity  of  the  electron  along  the  energy  axis  is 
defined  as 
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The  diffusion  coefficient  in  energy  space  is  written  analogously  to  that  in  geometric  space: 
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The  ratio  of  Ae/Ai  can  be  expressed  as  eEv,  since  the  electric  field  is  the  only  force  on  the 
electron  in  the  consideration  of  the  field  driven  flux.  Thus,  Eq  (A-3)  can  be  expressed  as 


D  =  i(veE)2A/ 


(A-4) 
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The  assumption  is  made  that  every  time  the  electron  collides  with  a  heavy  particle,  it  will  give 
up  all  of  its  excess  energy.  Doing  so  allows  At  to  be  expressed  as  1  Nm,  where  \m  is  the  momen¬ 
tum  transfer  collision  frequency.  Eq  (A-4)  can  then  be  written  as 
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2(eEfe 
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It  wiil  prove  to  be  advantageous  to  express  Eq  (A-5)  in  terms  of  E/N: 
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The  "velocity"  u  can  be  expressed  in  terms  of  the  diffusion  coefficient  as  (19): 
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Finally,  then,  the  current  density  in  energy  space  due  to  the  applied  field  can  be  expressed  as: 
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which  is  Eq  (5).  So  now  it  becomes  necessary  to  finite  difference  this  equation  in  energy  space. 
A  simple  backward  difference  expression  for  this  differencing  will  be  used: 
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It  will  be  useful  to  define  the  variable  r\k  such  that  it  represents  the  number  density  of  elec¬ 


trons  per  unit  energy  in  bin  k: 
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Doing  so  will  allow  the  following  expressions  to  be  redefined  as 


(  « 
nk  + 1 

y.  \ 

nk 

A £* 

^  Ae^  +  j 

_ 

A tkj 

2 

de 

ii 

si 

*- 

+ 

i 

si 

(A- 10a) 

(A- 10b) 


The  motivation  for  doing  this  lies  in  the  fact  that  the  electron  number  density  in  each  energy 
bin  needs  to  be  weighted  by  the  width  of  that  bin.  Using  Eq  (A- 10a)  and  (A- 10b),  J/fc)  can  be 
expressed  as 
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In  a  like  manner,  J ){k  -  1 )  can  be  written  as 
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Substituting  Eq  (A-ll)  and  (A-12)  into  Eq  (A-8),  and  grouping  into  like  terms  of  Tl*,  tj*  +  1. 


and  rji.,,  the  coefficients  can  be  written  as 
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where  the  prime  designates  that  these  are  the  coefficients  for  the  rit+1  and  r\k  terms  respectively. 
Flux  due  to  Elastics 

Rockwood  gives  the  electron  current  density  in  energy  space  due  to  elastic  energy  transfer  as 
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The  finite  differenced  expression  of  Eq  (A-14)  as  given  by  Rockwood  (14:2349)  allowed  A*’s 


that  could  be  negative,  thus  causing  numerical  instabilities  in  the  computation  of  the  solution. 
Following  the  guidelines  of  Greene  (7:2948),  the  following  approach  is  used: 
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where  to  represent  -dJjde  the  first  term  of  Eq  (A-15)  is  backward  differenced,  the  second  term 

is  forward  differenced,  and  the  last  is  backward  differenced.  Thus,  the  following  expressions  are 
used: 
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Again,  grouping  into  like  terms  of  tj*,  T]i+1,  and  T|4_„  the  coefficients  can  be  written  as 
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where  the  prime  designates  that  these  are  the  coefficients  for  the  T|i+1  and  terms  respectively. 

Putting  it  all  together 

Moving  back  to  the  original  representation  for  nk  and  nk  +  u  the  coefficients  become 
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In  practice,  it  has  been  found  that  the  use  of  these  equations  for  the  case  of  a  zeroed  electric 
field  results  in  the  distribution  obtaining  the  Maxwellian  steady  state  solution  as  expected.  In 
these  cases,  it  has  been  verified  that  the  average  energy  of  the  distribution  is  3/2  the  gas  tempera¬ 
ture.  However,  when  Eq  (17a)  adn  (17b)  are  used  in  conjunction  with  a  non-zero  electric  field, 
the  Bk  coefficients,  which  govern  the  rate  of  demotion  from  one  bin  to  the  next  lower  bin,  are  too 
large,  and  result  in  a  distribution  which  falls  off  much  faster  than  it  should.  Therefore,  in  the 
present  version  of  Megaboltz,  the  Rockwood  recoil  expressions  rederivcd  for  the  non-uniform 
energy  axis  are  used  when  an  electric  field  is  present,  while  Eq  (17)  is  used  with  no  electric  field 
present.  The  Rockwood  recoil  term  is  Eq  (15a)  expressed  as 
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Again  using  Eq  (A- 10a)  and  (A- 10b),  this  equation  can  be  finite  differenced  according  to  Eq 
(A-9).  Collecting  the  coefficients  of  like  electron  number  density  («*_,,  nk,  nk  +  l),  the  terms  for 
ak_  i,  (ak  +  bk),  and  6t  +  1  are  established.  The  resulting  expressions  are  given  as 
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Appendix  B 

Derivation  of  electron-electron  collision  term  for  a  non-uniform  energy  axis 


Let  an  energy  axis  be  defined  such  that  the  energy  of  bin  k,  et,  is  defined  at  the  center  of  the 
bin  with  a  width  Aet.  Using  this  definition,  nk  represents  the  number  density  of  electrons  with 

Ac,  Ac, 

energy  between  et  —  —  and  ek  +  — .  The  time  rate  of  change  of  the  electron  number  density  at  et 
is  given  by 
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Now,  let  the  electron  number  density  at  the  boundary  between  two  bins  to  defined  as 
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Similarly,  the  energy  derivative  of  the  electron  number  density  at  the  boundary  is  defined  as 
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Eq  (7)  is  finite  differenced  at  the  boundary  of  the  bins  k  and  k+1,  thus  becoming 
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A  similar  differencing  at  the  boundary  of  bins  k  and  k-1  results  in 


1k  - 1/2 


=  CX{ 


fp.+p.-l'l 

1 

( nk-nk_x  \ 

l  2  J 

_  2 £*  _  1/2 
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J? 

i 

1 

^Q*  +  Q*-i  Y  ^k-ink  +  Asknk_^ 


A 


Ae*_i  +  Ae* 


(B-5) 


Replacing  the  integrals  contained  in  Eq  (7)  with  summations. 


Pt  —  Elnl&£lHkt 

+  e,E7wn,A£,(l-H,,)]  (B-6) 

0*  =  (B-7) 

where  Hm  =  {° 

Substituting  Eq  (B-6)  and  (B-7)  into  (B-5),  and  collecting  the  coefficient  of  nk_l,ai_1  as  used 
in  Eq  (B-l)  is  determined.  This  coefficient  becomes 


ak- i  ~  ^  Ac  {fe  + 


3 

‘  Ae*  Y 

Z‘Uk~'  2 

t  Ae*-i  +  A£* 
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(B-8) 
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I 
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Defining  a\_ j  such  that 


s 

a  it-i  =  XA*.,/, 

the  coefficient  At ,  is  written  as 


(B-9) 


A,,  =  o^{(C.H,»u+^X.) 


e'“‘'2 


^£*+1 


Yl 


yj 


A£*  +  A£*  +  1_ 

+  [e*  +  id  -fW  +  e4(l  -Ht>/)](e7,/2Mi)}  (B-10) 

This  coefficient  for  nk  can  be  interpreted  as  the  rate  at  which  electrons  with  energy  ek  are 


promoted  to  energy  Ei+1.  This  promotion  is  accomplished  by  an  electron  with  energy  £,  being 
demoted  to  energy  e(_j.  An  expression  similar  to  Eq  (B-10)  can  be  derived  as  the  coefficient  of 
ni  +  l  by  using  Eq  (B-4).  When  transformed  into  the  k*  bin,  this  coefficient  becomes  BM,  and  is 
interpreted  as  the  rate  at  which  electrons  with  energy  ek  are  demoted  to  energy  e*  _ ,  by  electrons 
with  energy  e,  going  to  £/+i. 

Boundary  conditions  are  applied  to  Eq  (B-10)  and  the  corresponding  equation  for  Bk  l  so  that 


conservation  of  particles  is  strictly  observed.  These  boundary  conditions  given  by  Eq  (38),  and 
physically  mean  that  no  electron  can  be  demoted  from  the  lowest  energy  bin  when  k  =  1  ,  or 
promoted  from  the  highest  energy  bin  when  k  =  S  .  When  Eq  (B-10)  and  the  corresponding 
equation  for  Bt ,  are  programmed,  energy  is  not  conserved.  Additionally,  these  coefficients  will 
not  drive  the  distribution  toward  a  Maxwellian  over  time  as  they  should  when  only  electron  - 
electron  collisions  are  considered.  In  order  for  the  solution  method  to  strictly  observe  these  last 
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two  properties,  the  A*  ,’s  and  B*  /s  must  be  properly  constructed. 
Conservation  of  energy  can  be  expressed  as 


1  1 
dt  • 


oo 

Je/i(e)de 


t(e)de 


S  d/i. 

I 

*  =  i  * 

s 

I  "kte-k 

*=i 


(B-l  1) 


Using  Eq  (B-l)  as  the  time  derivative  of  the  electron  number  density,  and  the  definitions  of 
the  coefficients  a\.u  a'k,  b\,  and  b'k+l,  as  given  by  Eq  (37),  Eq  (B-l  1)  becomes 


Xrt*W/[(£*  +  i  ,/  ~  B 

k.l 


(B-12) 


This  condition  will  be  ensured  by 


b*,  = 

V  “  fc*  - 1  / 


(B-l  3) 


In  order  to  determine  the  constraint  on  the  matrices  caused  by  the  need  to  obtain  the  Maxwell¬ 
ian  in  the  steady  state,  the  microscopic  interpretations  of  Alt,  and  Bkl  are  used.  Once  the  steady 
state  Maxwellian  has  been  achieved,  then  the  flux  from  bin  k  upward  to  bin  k+1  must  be 
matched  exacdy  by  the  downward  flux  from  bin  k+1  to  bin  k.  Thus, 


f  Ae,  )  f  Ae;_, 

\,inkni  =  ^k  +  i,t-\nk  +  int-i  ~7Z 

V  ^k  y  V  Gt'k  + 1 


(B-14) 


Using  Eq  (B-13)  to  express  Bt  +  I in  terms  of  the  A  matrix,  Eq  (B-14)  becomes 
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(B-15) 


Since  the  Maxwellian  distribution  defines  the  electron  number  density  nk  as 


nk  =  Ce,/2exp 


(B-16) 


the  two  number  density  ratios  can  be  rewritten  and  the  final  expression  for  Eq  (B-14)  becomes 


A-i,*  +  : 


r  f -e  v 

t/  t-i  _  i 


£*  + 1  —  £*  y 


V 


£*  +  i 
£* 


y/2/  y/2 

fcy-1 

V  £'  J 


x 


exri 


(A£i  +  A£;-1)-(A£,4l-A£,) 
2T 


(B-17) 


In  practice,  Eq  (B-10)  is  used  to  define  the  part  of  the  A  matrix  on  and  below  the  diagonal,  as 
well  as  the  elements  just  above  the  diagonal.  Eq  (B-17)  is  then  used  to  define  the  rest  of  the  A 
matrix.  The  first  column  and  last  row  of  the  A  matrix  is  set  to  zero,  which  will  ensure  the  con¬ 
servation  of  particles  as  demanded  by  Eq  (38),  then  the  entire  B  matrix  is  constructed  by  the  use 
of  Eq  (B-13). 
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The  current  version  of  MEGABOLTZ  is  designed  to  run  under  UNIX  version  4.3,  with  the 
IMSL  v9.2  library  and  METALIB  library  installed.  There  are  currently  two  calls  to  the  IMSL 
library  for  the  L-U  decomposition  and  substitution,  while  all  graphing  functions  are  accom¬ 
plished  by  the  METALIB  library.  In  order  to  run,  MEGABOLTZ  expects  some  files  to  be  in  the 
same  directory  as  the  executable  code.  These  files  consist  of  the  following: 

EILE  DESCRIPTION 

megaboltz  executable  Boltzmann  solver  code 

input.com  ASCII  data  file  describing  the  physical  conditions  under  the 

solution  is  to  be  obtained 

*.crs  ASCII  cross  section  data  file,  ie.,  ar.crs  for  Argon  or  h2.crs  for 

Hydrogen.  This  file  is  explained  in  detail  in  Appendix  D. 

The  input  file,  input.com,  is  shown  in  Figure  36.  The  E/N  select  switch  allows  the  user  to 
chose  if  the  calculation  is  for  one  E/N  or  for  a  list  of  E/N’s.  If  a  list  of  E/N’s  is  used,  the  list  may 
be  either  equally  spaced  (choice  3),  or  contained  in  an  external  file  of  E/N’s  called  "eovem.dat", 
(choice  2).  If  the  equally  spaced  option  is  used,  the  number  of  data  points  is  calculated  by  [(Start 
E/N  -  End  E/N)  /  delta  E/N  ].  For  example,  suppose  you  are  interested  in  running  2  to  30  Td,  by 
2  Td  increments.  Select  the  start  E/N  value  to  be  30.00d-17,  and  the  end  value  to  be  2.00d-17. 
The  number  of  datapoints  is  [(30-2)/2]  =  14.  Notice  that  the  actual  number  of  datapoints  will  be 
one  more  than  this,  or  15.  The  start  value  is  given  as  30  Td  instead  of  2  Td  because  MEGA¬ 
BOLTZ  expects  the  higher  value  to  be  given  first.  This  choice  was  driven  by  the  quicker 
convergence  at  the  higher  E/N  values.  The  initial  guess  distribution  for  the  list  of  E/N’s  will  be  a 
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Maxwellian  at  the  start  E/N  value.  For  each  next  E/N  value,  the  initial  guess  will  be  the 
converged  solution  for  the  previous  E/N.  Doing  things  in  this  order  speeds  the  calculation  of  the 
distribution  at  each  E/N. 


* —  Input  Parameters 
E/N  Select: 

1 - Single  Point  ->  Solve  for  START  E/N  value  only 

2- Read  Ext  File  ->  Solve  for  E/N’s  in  File  ’eovem.dat’ 

3- Iterate  ->  Solve  for  equally  spaced  E/N’s  below 

START  E/N  (V-CmA2):  1.0D-17 

END  E/N  :  40.00D-17 
DATA  Points  :  8 


E/N  Choice  (1,2,3) :  1 

Gas  Temperature  (K):  3.d2 
Gas  Pressure  (Torr):  7.6d2 
e-  Num  Density  (CmA-3):  1.0D1 3 


* —  Mesh  Parameters 

Time  Step  (S):  l.d-9 
Number  Of  Bins:  150 
Bin  Width  Select: 

1 - VARIABLE  ->  Variable  Bin  Width 

2- LOG  ->  Logarithmic  Bin  Width 

3- L  T'fEAR  ->  Constant  Bin  Width 
Max  Energy  (eV):  10.0 


Bin  Width  (1,2,3):  3 

* —  Program  Switches 
E-E  Collisions?  (Y/N):  N 
e-  Beam?  (Y/N):  N 
Source  current:  l.d-3 
Source  location:  30.0 
Save  Data  As:  testdat 
Save  Graph  As:  test.plt 

*—  Gasmix  data  by  type  (5  maximum),  percent  (enter  as  XXX),  and  whether 
you  want  to  view  its  cross-sections  (enter  either  Y  or  N) 

_Gas _ Percent _ View?_ 

ar  99.9  N 

f  0.1  N 

EOF  0.0  N 

Figure  36.  Sample  Input  File  "input.com"  for  MEGABOLTZ 
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The  mesh  parameters  determine  the  energy  grid  along  which  the  solution  will  be  determined. 
For  a  uniform  grid,  the  number  of  energy  bins  and  the  maximum  energy  will  determine  the 
energy  bin  width.  In  the  case  of  a  Log  grid,  this  maximum  energy  has  no  meaning.  Instead,  the 
program  will  prompt  the  user  for  an  alpha  value  (Eq  (64)).  The  maximum  energy  is  determined 
by  the  sum  of  each  energy  bin  defined  by  this  parameter.  By  varying  alpha,  the  maximum 
energy  can  be  selected.  The  code  will  loop  through  this  process  until  the  user  agrees  to  procede 
with  the  resulting  maximum  energy.  The  Variable  grid  allocation  method  is  presently  under 
development,  and  its  use  is  presently  not  recommended.  The  time  scale  to  be  used  is  also  set  in 
this  block.  This  sets  the  time  step  of  integration  for  the  coupled  set  of  differential  equations, 
given  in  Eq  (1 1)  and  Eq  (23).  This  value  should  be  at  least  as  small  as  the  characteristic  time 
scale  for  electron  interaction,  given  by  1/v,  where  v  is  the  collisional  frequency  for  the  dominate 
process. 

The  file  names  under  which  the  data  and  graphs  should  be  saved  are  listed  in  the  next  block, 
along  with  the  swithes  governing  electron-electron  collisions  and  electron  beams.  If  a  duplicate 
file  name  is  used,  MEGABOLTZ  will  prompt  the  user  for  a  new  file  name.  The  final  block  lists 
the  gases  to  include  in  the  calculation.  The  gases  listed  in  this  section  must  be  spelled  exactly  as 
the  prefix  to  the  .crs  extension  for  the  corresponding  gas.  For  example,  if  the  cross  section  file  is 
ar.crs,  then  the  gas  listed  here  is  ar,  not  Ar.  After  the  final  gas,  EOF  must  be  placed  as  the  end  of 
file  marker,  along  with  a  gas  concentration  of  0.0%.  This  an  important  step,  because  without  it 
MEGABOLTZ  may  be  able  to  compute  a  converged  solution  that  is  incorrect. 

Once  the  solution  has  converged,  MEGABOLTZ  will  save  all  the  transport  and  energy  bal¬ 
ance  data  to  the  file  specified  in  "input.com".  The  solution  will  be  saved  in  a  file  named  "N.dat”. 
The  graph  of  this  file  will  be  saved  to  the  file  specified  in  "input.com".  A  listing  of  all  files 
written  by  MEGABOLTZ  is  shown  below. 

EILE  DESCRIPTION 
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testdat 

testplt 


N.dat 

xsmod.dat 


Name  of  ASCII  file  specified  by  user  in  "input.com"  containing 
transport  and  energy  balance  data  for  each  E/N  calculation. 

Name  of  file  specified  by  user  in  "input.com"  containing  the 
plot  file  of  the  distribution.  This  file  is  constructed  so  as  to  print 
out  on  a  laser  printer  when  used  with  the  command  mit  "test.plt" 
>  lpr  -Pimagen,  from  the  Unix  sytem  prompt. 

Name  of  ASCII  file  containing  the  electron  energy  distribution 
of  the  converged  solution  for  each  energy  bin. 

Name  of  ASCII  file  containing  cross  section  information.  This 
file  is  written  as  each  data  point  from  the  cross  section  file  *.crs 
is  read.  It  is  used  mainly  for  debugging  *.crs  files.  (Appendex 
D) 


The  sample  output  file  "test.dat"  is  shown  below: 


Megaboltz  Output  Data 
Bins  Binwidth  Timestep 
150  0.1000  1.00E-09 
Nl/N  :  0.00E+00 
Ne  final:  8.75437E+12 

E/N  Rion  Rion2  Rexc  Rsuper  Rattach 

1.0000E-16  1.2312E-I7  O.OOOOE+OO  5.9674E-12  O.OOOOE+OO  O.OOOOE+OO 

E/N  <E>  Echar  Vd  Df  E  Balance 

1.000D-16  5.371D+00  7.819D+00  9.221D+05  2.946D+03  3.784D-01 

Megaboltz  Energy  Calculations 
Energy  Loss  Mechanisms 

E/N  Sys  Loss  El  Losses  Exc  Losses  Ion  Losses  Atch  Losses 

1.0000E-16  2.2655E+09  6.2973E+08  1.7045E+09  3.259E+04  0.0000+00 

Energy  Gain  Mechanisms 
E/N  Sys  Gains  J.E  Super  e  Beam 

1.0000E-16  2.2569E+09  2.2569E+09  O.OOOOE+OO  O.OOOOE+OO 

Figure  37.  Sample  Output  File  "test.dat"  From  MEGABOLTZ 


Finally,  it  must  be  stressed  that  MEGABOLTZ  expects  the  parameters  listed  in  "input.com" 
to  be  in  the  proper  place.  If  the  template  given  for  "input.com"  is  followed  exactly,  no  problems 
will  occur.  However,  MEGABOLTZ  has  not  been  "idiot-proofed",  and  it  is  adviseable  to  keep  a 
backed  up  copy  of  this  file.  Even  experienced  users  have  found  it  to  be  advantageous. 
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Just  as  MEGABOLTZ  expects  the  "input.com"  file  to  properly  formatted,  it  expects  each 
cross  section  file  to  be  properly  formatted.  The  purpose  of  this  appendex  is  to  familiarize  the 
user  with  the  construction  of  the  cross  section  file  to  the  extent  that  modification  can  be  com¬ 
pleted  relatively  painlessly.  Again,  having  a  backed  up  copy  of  each  cross  section  file  is  an 
extremely  good  idea.  Figure  38  shows  the  cross  section  data  file  for  Argon. 

The  first  line  gives  just  the  gas  name.  This  title  is  not  read  by  MEGABOLTZ,  so  it  can  be  in 
any  format.  The  next  line  contains  three  pieces  of  information.  The  first  four  spaces  are  dedi¬ 
cated  to  the  molecular  weight  for  the  gas,  while  the  next  three  are  used  to  designate  how  many 
cross  sections  are  included  in  the  file.  This  must  be  an  integer  between  1  and  999  (which  would 
be  a  lot  of  cross  sections!).  The  last  number  is  the  ionization  threshold  in  eV  for  the  gas,  how¬ 
ever  it  is  not  read  by  MEGABOLTZ.  Each  line  is  ended  by  a  forward  slash  (/)  that  must  be  out¬ 
side  of  the  defined  block  area. 

Each  cross  section  table  within  the  cross  section  file  consists  of  three  "lines",  except  the  first 
table  which  is  the  momentum  transfer  table  and  consists  on  only  two  "lines".  These  lines  are  not 
to  be  confused  with  lines  of  ASCII,  but  rather  have  no  particular  maximum  length.  The  end  of 
one  of  these  line  is  defined  by  a  forward  slash  (/).  The  first  line  of  each  table  begins  with  a  three 
space  block  dedicated  to  the  number  of  data  point  pairs  (energy,  cross  section)  in  the  table.  The 
number  occupying  this  block  must  be  an  integer.  The  next  12  spaces  are  in  a  reserved  block 
which  designates  the  type  of  process  the  table  contains.  For  example,  this  identifier  could  be 
"elastic",  "inelastic",  "ionization",  or  "attachment",  where  the  quote  marks  are  not  used.  In  the 
process  of  running,  MEGABOLTZ  will  look  for  these  identifying  words  as  labels,  thus  they  must 
be  spelled  correctly,  with  no  capitals.  The  next  15  spaces  are  dedicated  to  a  block  that  describes 
the  origin  of  the  data.  The  quotes  are  included  as  part  of  the  15  total  spaces  of  the  block.  The 
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next  15  spaces  are  dedicated  to  another  block,  which  can  be  used  as  a  continuation  of  the  pre¬ 
vious  one.  Again,  the  quotes  arc  included  as  part  of  the  15  space  total.  The  final  variable  on  this 
line  represents  the  scale  of  the  data  points.  This  occupies  4  spaces  (XXX.X).  This  is  used  as  a 
multiplier  for  the  cross  sectional  values,  and  is  usually  1 .0.  This  line  is  ended  with  a  forward 
slash  (/)  that  is  outside  of  the  defined  block  area,  however,  remember  that  Fortran  77  can  only 
read  72  characters  per  line. 


AR 

40.5  15.759/ 

54  elastic  ’FAUST.J.PHYS,”  1977,30,6 1-72’  1.0/ 

0.000,0.000  0.014,3.880  0.017,3.560  0.020,3.280  0.025,2.890 
0.030,2.570  0.035,2.290  0.040,2.050  0.050,1.662  0.060,1.357 
0.070,1.114  0.080,0.916  0.090,0.754  0.100,0.621  0.110,0.511 
0.120,0.420  0.130,0.348  0.140,0.284  0.150,0.233  0.170,0.161 
0.180,0.135  0.190,0.115  0.200,0.101  0.210,0.092  0.220,0.086 
0.230,0.085  0.240,0.087  0.050,0.091  0.260,0.098  0.280,0.120 
0.300,0.151  0.320,0.188  0.325,0.206  0.400,0.317  0.500,0.504 
0.650,0.792  0.800,1.050  1.000,1.370  1.200,1.660  1.500,2.050 
1.700,2.330  2.000,2.700  2.500,3.430  3.000,4.200  4.000,5.700 
5.000,7.600  7.600,12.50  10.00,16.61  15.00,16.14  20.00,12.50 
40.00,6.0  60.00,4.6  80.00,4.1  100.00,4.0  / 

27  inelastic  ’AR  El’  ’  ’  1.0  / 


0.  , 

0.  11.6  ,  0. 

11.65  , 

.0037 

11.7 

,  .0074  11.75  , 

.0104  11.8 

,  .0134 

11.85 

,  .0141  11.9  , 

.0149  12. 

,  .0163 

12.1 

,  .0168  12.2  , 

.017  12.4 

,  .018 

12.6 

,  .0188  12.8  , 

.0196  13. 

,  .0205 

13.2 

,  .0213  13.4  , 

.022  14. 

,  .0225 

20. 

.  .0271  30.  , 

.0275  40. 

,  .0275 

50. 

,  .0233  60.  , 

.0203  70. 

,  .0181 

80. 

,  .0163  90.  , 

.0149  100. 

,  .0137 

11.6,  1.0,  0.0  / 


23  inelastic  ’AR  E2’  ’  ’  1.0  / 

0.  ,  0.  11.8  ,  0.  11.85  ,  .0022 


11.9 

,  .0078  12.  , 

.0134 

12.1 

,  .0172 

12.2 

,  .0214  12.4  , 

.0295 

12.6 

,  .0411 

12.8 

,  .0516  13.  , 

.0643 

13.2 

,  .0753 

13.4 

,  .0778  13.6  , 

.078 

14. 

,  .0784 

30. 

,  .0944  40.  , 

.0912 

50. 

,  .0774 

60. 

,  .0675  70.  , 

.0600 

80. 

,  .0542 

90. 

,  .0494  100.  , 

.0455 

/ 

11.8, 1.0, 0.0/ 

24  inelastic  ’AR  E3’ 

» 

* 

1.0  / 

0. 

,  0.  13.2  ,  0 

.  13.4  ,  . 

0013 

13.5 

,  .0018  13.6  , 

.0065 

l  13.7 

,  .0133 

13.8 

,  .0206  14.  , 

.0359 

15. 

,  .1323 

16. 

,  .2826  17.  , 

.465 

18.  , 

,  .645 

106 


19.  ,  .82  20.  ,  1.16  22.  ,  1.27 

26.  ,  1.48  30.  ,  1.56  40.  ,  1.48 

50.  ,  1.256  60.  ,  1.095  70.  ,  .9740 

80.  ,  .8790  90.  ,  .8024  100.  ,  .7390  / 
13.2, 1.0, 0.0/ 

38 ion  ’ART  ’KIEFFER,78’  1.0  / 
0.  ,  0.  15.759 ,  0.  16.  ,  .02023 
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21. 
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23.  , 
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24. 
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25.  , 
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26.  , 
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28. 
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30.  , 

1.803 

32.  , 
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34. 

,  2.111 
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2.243 
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2.331 

40. 

,  2.393 

42.5  , 

2.446 

45.  , 

2.49 

50. 

2.534 

55.  , 

2.595 

60.  , 

2.657 

65. 

,  2.727 

70.  , 

2.771 

75.  , 

2.815 

80. 

,  2.841 

85.  , 

2.85 

90.  , 

2.859 

100. 

,  2.85 
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2.824 

130.  , 

2.762 

145. 

,  2.709 

160. 
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180.  , 

2.516 

200. 

,  2.393 

/ 

15.759 

,  1.0/ 

--EOR- 

* 

--EOF- 

- 

Figure  38.  Cross  Section  Data  File  for  Argon 


The  second  line  of  each  table  actually  consists  of  many  lines  of  numbers,  representing  the 
pointwise  description  of  the  cross  section.  These  data  points  are  formatted  as  (Energy,  Cross 
section),  where  the  energy  is  in  units  of  eV,  and  the  cross  section  default  is  units  of  10E-16 
cmA2.  This  default  can  be  modified  with  the  scale  variable  defined  earlier.  There  is  a  comma 
between  each  number  of  the  data  pair,  while  at  least  one  space  separates  each  data  pair  from  each 
other.  Once  each  of  the  data  pairs  in  each  table  is  listed,  a  forward  slash  (/)  designates  the  end  of 
the  line. 

The  third  line  is  included  only  for  cross  section  tables  after  the  momentum  transfer  table.  The 
format  of  this  third  tine  depends  on  whether  the  table  is  for  an  "inelastic"  process  or  not.  If  the 
table  is  for  "inelastic",  then  the  third  tine  consists  of  three  blocks.  If  the  table  is  for  "ionization" 
or  "attachment",  then  this  third  line  is  constructed  of  only  two  blocks.  Each  of  these  blocks  has 
no  specific  length,  rather,  they  are  separated  by  commas.  The  first  block  of  both  types  of  tines 
designates  the  threshold  energy  for  the  process,  in  eV.  The  second  number  in  the  third  line  rep- 
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resents  the  ratio  of  the  number  of  gas  particles  to  be  affected  by  the  process  to  the  number  of  gas 
particles  in  the  ground  state.  Thus,  this  number  must  be  between  0.0  and  1.0.  When  this  number 
is  multiplied  by  the  total  number  of  gas  molecules,  the  result  is  the  number  of  molecules  affected 
by  the  state.  For  any  process  affecting  ground  state  molecules,  this  number  will  be  1.0.  This 
includes  ground  state  ionization  and  excitation,  as  well  as  a  process  such  as  attachment.  By  mak¬ 
ing  the  value  of  this  variable  0.0  for  any  of  these  processes,  that  particular  process  can  be 
"tumed-off '  in  this  manner.  For  "inelastic"  processes,  the  third  number  designates  the  number  of 
superelastic  molecules  to  consider.  This  number  is  a  ratio  of  the  number  of  molecules  in  the 
upper  excited  state  of  the  excitation  process,  to  the  number  in  the  lower  state  of  the  same  pro¬ 
cess.  Consider  an  excitation  process  from  the  ground  state  to  some  upper  excited  state  j.  If  only 
1  out  of  every  10,000  molecules  is  in  the  state  j,  then  this  number  would  be  1.0E  5.  In 
calculating  the  effects  of  superelastic  collisions,  the  code  will  use  this  number  in  determining  the 
number  of  collisional  partners  for  the  electrons.  This  line  is  ended  with  a  forward  slash  (/). 

The  selection  of  these  parameters  formatted  in  this  manner  may  seem  to  be  ackward,  but  after 
practice  their  use  will  become  second  nature  and  comfortable.  The  hard  structure  of  the  cross 
section  files  is  very  strict,  and  very  unforgiving.  In  order  to  help  the  user  with  the  formatting  of 
the  files,  MEGABOLTZ  will  recreate  the  cross  section  file  as  it  reads  it.  This  file,  called 
"xsmod.dat",  is  created  each  time  MEGABOLTZ  runs,  and  does  not  have  to  be  in  existence  prior 
to  running  the  code.  This  file  actually  contains  two  files:  the  first  file  is  the  cross  section  file 
rewritten,  while  the  second  file  consists  of  the  interpolated  cross  sections  for  up  to  five  tables. 
Many  times  MEGABOLTZ  errors  can  be  traced  to  incorrectly  formatted  cross  section  files, 
"xsmod.dat"  is  designed  to  help  iron  out  these  errors. 
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